home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-19 / trjmkr.zip / MATH20.DOC < prev    next >
Text File  |  1992-08-09  |  109KB  |  3,607 lines

  1.  
  2.  
  3.  
  4.  
  5.  
  6.  
  7.  
  8.  
  9.  
  10.                         TRAJECTORY MAKER
  11.  
  12.                         A Trajectory Simulation Tool
  13.                         ────────────────────────────
  14.  
  15.  
  16.  
  17.  
  18.  
  19.  
  20.  
  21.  
  22.  
  23.  
  24.  
  25.  
  26.  
  27.  
  28.  
  29.  
  30.  
  31.  
  32.  
  33.  
  34.  
  35.  
  36.  
  37.                         
  38.                                     ALGORITHM SUPPLEMENT
  39.  
  40.                         
  41.  
  42.  
  43.  
  44.  
  45.  
  46.  
  47.  
  48.  
  49.  
  50.                         
  51.  
  52.  
  53.  
  54.                                                              ORBIT HORIZONS
  55.                                                                94 Promenade
  56.                                                    Irvine, California 92715
  57.                                                                      U.S.A.
  58.  
  59.  
  60.  
  61.  
  62.                                 COPYRIGHT NOTICE
  63.  
  64.      This document and software package  consisting of the Trajectory Maker
  65.      and Trajectory Scape computer programs  are copyrighted (C) 1991, 1992
  66.      by H.B. Reynolds, President, ORBIT HORIZONS.  All rights are reserved.
  67.      No  part   of  this   publication  may   be  reproduced,  transmitted,
  68.      transcribed, stored  in any retrieval  system, or  translated into any
  69.      language by any means without  the express written permission of ORBIT
  70.      HORIZONS, 94 Promenade, Irvine California 92715, USA.
  71.  
  72.  
  73.  
  74.  
  75.  
  76.  
  77.  
  78.  
  79.  
  80.  
  81.  
  82.  
  83.  
  84.  
  85.  
  86.  
  87.  
  88.  
  89.  
  90.  
  91.  
  92.  
  93.  
  94.  
  95.  
  96.  
  97.                           FIRST EDITION/FIRST PRINTING
  98.  
  99.                                     July 1992
  100.  
  101.  
  102.  
  103.  
  104.  
  105.  
  106.  
  107.  
  108.                                                                          ii
  109.  
  110.  
  111.  
  112.  
  113.                                 LIMITED WARRANTY
  114.  
  115.      THIS  ALGORITHM  SUPPLEMENT,  THE  USER'S  GUIDE,  AND  THE  PROGRAMS,
  116.      "TRAJECTORY MAKER" AND  "TRAJECTORY SCAPE" ARE  SOLD "AS IS",  WITHOUT
  117.      WARRANTY  AS  TO  THEIR PERFORMANCE,  OR  FITNESS  FOR  ANY PARTICULAR
  118.      PURPOSE.  THE  ENTIRE RISK  AS TO THE  RESULTS AND  PERFORMANCE OF THE
  119.      PROGRAMS IS ASSUMED BY THE USER.
  120.  
  121.      HOWEVER, ORBIT HORIZONS WARRANTS THE MAGNETIC DISKETTE(S) ON WHICH THE
  122.      PROGRAM IS RECORDED  TO BE FREE  FROM DEFECTS IN  MATERIALS AND FAULTY
  123.      WORKMANSHIP UNDER NORMAL USE FOR A PERIOD OF NINETY DAYS FROM THE DATE
  124.      OF PURCHASE.   IF DURING  THIS NINETY  DAY PERIOD  THE DISKETTE SHOULD
  125.      BECOME  DEFECTIVE,  IT  MAY  BE  RETURNED  TO  ORBIT  HORIZONS  FOR  A
  126.      REPLACEMENT WITHOUT CHARGE, PROVIDED YOU SEND PROOF OF PURCHASE OF THE
  127.      PROGRAM.
  128.  
  129.      YOUR SOLE AND EXCLUSIVE  REMEDY IN THE EVENT  OF A DEFECT IS EXPRESSLY
  130.      LIMITED TO REPLACEMENT OF THE  DISKETTE AS PROVIDED ABOVE.  IF FAILURE
  131.      OF  A DISKETTE  HAS RESULTED  FROM ACCIDENT  OR ABUSE,  ORBIT HORIZONS
  132.      SHALL HAVE NO  RESPONSIBILITY TO REPLACE THE  DISKETTE UNDER THE TERMS
  133.      OF THIS LIMITED WARRANTY.
  134.  
  135.      IN NO  EVENT SHALL ORBIT  HORIZONS BE  LIABLE TO YOU  FOR ANY DAMAGES,
  136.      INCLUDING BUT NOT LIMITED TO ANY  LOST PROFITS, LOST SAVINGS, OR OTHER
  137.      CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THESE
  138.      PROGRAMS AND  ACCOMPANYING MATERIALS EVEN  IF ORBIT  HORIZONS HAS BEEN
  139.      ADVISED OF THE  POSSIBILITY OF SUCH  DAMAGES, OR FOR  ANY CLAIM BY ANY
  140.      OTHER PARTY.
  141.  
  142.      SOME STATES DO NOT ALLOW THE  LIMITATION OR EXCLUSION OF LIABILITY FOR
  143.      INCIDENTAL  OR  CONSEQUENTIAL  DAMAGES,  SO  THE  ABOVE  LIMITATION OR
  144.      EXCLUSION  MAY NOT  APPLY TO  YOU.  THIS  WARRANTY GIVES  YOU SPECIFIC
  145.      LEGAL RIGHTS, AND YOU MAY ALSO HAVE OTHER RIGHTS WHICH VARY FROM STATE
  146.      TO STATE.
  147.  
  148.      THIS WARRANTY SHALL BE CONSTRUED,  INTERPRETED AND GOVERNED BY LAWS OF
  149.      THE STATE OF CALIFORNIA.
  150.  
  151.      YOU AGREE TO THESE TERMS BY YOUR DECISION TO USE THIS SOFTWARE.
  152.  
  153.  
  154.  
  155.  
  156.  
  157.  
  158.  
  159.  
  160.  
  161.  
  162.                                                                               iii
  163.  
  164.  
  165.  
  166.  
  167.  
  168.  
  169.                                                                        iii
  170.  
  171.  
  172.  
  173.  
  174.  
  175.  
  176.  
  177.                                    TABLE OF CONTENTS
  178.  
  179.  
  180.                1  INTRODUCTION ...................................... 1-1
  181.  
  182.                2  EQUATIONS OF MOTION ............................... 2-1
  183.  
  184.                3  GRAVITY PERTURBATION MODEL ........................ 3-1
  185.  
  186.                4  DRAG PERTURBATION MODEL ........................... 4-1
  187.                    4.1  Density Profile ............................. 4-4
  188.                    4.2  Wind Profile ............................... 4-14
  189.  
  190.                5  NUMERICAL INTEGRATION SCHEME ...................... 5-1
  191.  
  192.                6  COORDINATE SYSTEMS AND TRANSFORMATIONS ............ 6-1
  193.  
  194.                    6.1  Time Considerations ......................... 6-1
  195.                    6.2  Orbital Elements ............................ 6-2
  196.                    6.3  Geospherical Inertial Coordinates ........... 6-6
  197.                    6.4  Earth-Fixed-Geocentric Coordinates .......... 6-6
  198.                    6.5  Local Topocentric Coordinates ............... 6-7
  199.                    6.6  Common Radar Coordinates ................... 6-10
  200.                    6.7  Geodetic Coordinates ....................... 6-10
  201.  
  202.                7 BIBLIOGRAPHY ....................................... 7-1
  203.  
  204.  
  205.  
  206.                                     LIST OF TABLES
  207.  
  208.  
  209.                Table                                                Page
  210.  
  211.                 3-1  WGS-72 Earth Constants ........................ 3-1
  212.  
  213.                 3-2  Smithsonian Zonal Harmonic Coefficients ....... 3-2
  214.  
  215.                 4-1  Drag Coefficient vs. Mach Number .............. 4-3
  216.  
  217.                 4-2  Molecular Scale Temperature ................... 4-4
  218.  
  219.                 4-3  Atmospheric Density Profile ................... 4-6
  220.  
  221.                 5-1  Shanks 8-12 Numerical Integration Constants ... 5-3
  222.  
  223.  
  224.  
  225.  
  226.  
  227.  
  228.  
  229.                                                                        iv
  230.  
  231.  
  232.  
  233.  
  234.                                     LIST OF FIGURES
  235.  
  236.  
  237.                Figure                                                Page
  238.  
  239.                 3-1  Magnitude of Vector Differences Between the
  240.                      23rd and 2nd, 3rd, 4th and 5th Order Gravity
  241.                      Models ......................................... 3-7
  242.  
  243.                 3-2  Magnitude of Vector Differences Between the
  244.                      23rd and 6th, 7th, 8th and 9th Order Gravity
  245.                      Models ......................................... 3-8
  246.  
  247.                 3-3  Magnitude of Vector Differences Between the
  248.                      23rd and 10th, 11th, 12th and 13th Order
  249.                      Gravity Models ................................. 3-9
  250.  
  251.                 3-4  Magnitude of Vector Differences Between the
  252.                      23rd and 14th, 15th, 16th and 17th Order Gravity
  253.                      Models ........................................ 3-10
  254.  
  255.                 3-5  Magnitude of Vector Differences Between the
  256.                      23rd and 18th, 19th, 20th and 21st Order
  257.                      Gravity Models ................................ 3-11
  258.  
  259.                 3-6  Magnitude of Vector Differences Between the
  260.                      23rd and 19th, 20th, 21st and 22nd Order
  261.                      Gravity Models ................................ 3-12
  262.  
  263.                 3-7  Magnitude of Vector Differences Between the
  264.                      23rd and 2nd, 4th, 6th and 8th Order (Even)
  265.                      Gravity Models ................................ 3-13
  266.  
  267.                 3-8  Magnitude of Vector Differences Between the
  268.                      23rd and 3rd, 5th, 7th and 9th Order (Odd)
  269.                      Gravity Models ................................ 3-14
  270.  
  271.                 3-9  Magnitude of Vector Differences Between the
  272.                      23rd and 2nd, 6th, 10th and 14th Order
  273.                      Gravity Models ................................ 3-15
  274.  
  275.                 3-10 Magnitude of Vector Differences Between the
  276.                      23rd and 5th, 10th, 15th and 20th Order
  277.                      Gravity Models ................................ 3-16
  278.  
  279.                 4-1  1976 U. S. Standard Atmosphere Density
  280.                      Profile ........................................ 4-7
  281.  
  282.                 4-2  Abbreviated 1976 U. S. Standard Atmosphere
  283.                      Density Profile ................................ 4-8
  284.  
  285.  
  286.                                                                         v
  287.  
  288.  
  289.  
  290.  
  291.                                     LIST OF FIGURES
  292.  
  293.  
  294.                Figure                                                Page
  295.  
  296.                 4-3  Trajectory Maker Density Profile Model ......... 4-9
  297.  
  298.                 4-4  Cubic Spline Interpolation Error for 21 and
  299.                      48 Point Tables ............................... 4-10
  300.  
  301.                 4-5  Cubic Spline Interpolation Error for 21 and
  302.                      48 Point Tables ............................... 4-11
  303.  
  304.                 4-6  Cubic Spline and Linear Interpolation Error
  305.                      for 48 Point Table ............................ 4-12
  306.  
  307.                 4-7  Cubic Spline Interpolation Error for 48 Point
  308.                      Table ......................................... 4-13
  309.  
  310.                 4-8  East/North Wind velocity Components
  311.                      (Example 1) ................................... 4-15
  312.  
  313.                 4-9  East/North Wind velocity Components
  314.                      (Example 2) ................................... 4-16
  315.  
  316.                 5-1  Numerical Integration Error .................... 5-6
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  
  341.  
  342.  
  343.  
  344.                                                                        vi
  345.  
  346.  
  347.  
  348.  
  349.                Trajectory Maker                              Introduction
  350.                Algorithm Supplement                                  1-1
  351.                ──────────────────────────────────────────────────────────
  352.                
  353.  
  354.  
  355.                1.  INTRODUCTION
  356.  
  357.                The Trajectory Maker/Trajectory Scape software system is a
  358.                state of the art simulation and design tool for use on IBM
  359.                compatible personal computers.  Within minutes it predicts
  360.                precision trajectory ephemerides and displays their ground
  361.                traces on world map projections.
  362.  
  363.                The shapes  of various  orbit ground  traces can  be quite
  364.                bizarre, unexpected and difficult  to visualize because of
  365.                two motion  processes; namely  the Earth's  rotation about
  366.                its axis, coupled with the  motion of the orbiting object.
  367.                With  Trajectory  Maker,  visualization  of  this  complex
  368.                interaction is a breeze.
  369.  
  370.                But Trajectory Maker  gives much more  than this graphics'
  371.                visualization capability.  It  also creates high-precision
  372.                numerical  output  data,  in  several  formats,  for those
  373.                wanting more comprehensive analysis capabilities.
  374.  
  375.                Major  features of  the  software and  how  to use  it are
  376.                described in detail in the User's Guide.  The User's Guide
  377.                also  contains  some  historical  background  material,  a
  378.                glossary of  terms, definitions  and illustrations  to aid
  379.                the user with the external interface to the system.
  380.  
  381.                Transparent  to the  user is  an underlying  foundation of
  382.                sophisticated mathematical algorithms  and techniques that
  383.                are  employed  by  the  system  in  order  to  perform the
  384.                relevant  numerical  computations,   in  an  accurate  and
  385.                expedient manner.
  386.  
  387.                The user of  the system need  not necessarily be concerned
  388.                with  its  inner workings;  however,  an  understanding of
  389.                theory is sometimes helpful in understanding its practice.
  390.  
  391.                This paper presents a detailed description of mathematical
  392.                algorithms and  techniques used  by the  Trajectory Maker/
  393.                Trajectory Scape software system.
  394.  
  395.                In  addition to  presenting the  mathematical core  of the
  396.                software system, the  merits and limitations  of the algo-
  397.                rithms  are  discussed.  Illustrations  are  provided that
  398.                will aid in selecting various  aspects of the models, such
  399.                as truncation error, integration  step-size, etc., for the
  400.                particular application under consideration.
  401.  
  402.                The equations  of motion  that govern  the trajectories of
  403.  
  404.  
  405.  
  406.  
  407.                Trajectory Maker                              Introduction
  408.                Algorithm Supplement                                  1-2
  409.                ──────────────────────────────────────────────────────────
  410.                
  411.  
  412.  
  413.                Earth  orbiting  satellites  and  ballistic  missiles  are
  414.                summarized.
  415.  
  416.                The force models include perturbation effects due to Earth
  417.                oblateness, aerodynamic drag and atmospheric winds.  These
  418.                forces are developed and investigated.
  419.  
  420.                Explicit  closed-form  expressions  for  the gravitational
  421.                acceleration components  are derived  for an axi-symmetric
  422.                oblate  Earth  model.   The  model  contains  an arbitrary
  423.                number of zonal harmonic coefficients.  Effects of various
  424.                orders of the coefficients  on certain mission ephemerides
  425.                are investigated.
  426.  
  427.                The  aerodynamic drag  force  is defined,  the atmospheric
  428.                density profile is  modeled as a  function of altitude and
  429.                its relative  error is analyzed.   An innovative algorithm
  430.                for computing altitude is presented that is very fast, yet
  431.                produces highly accurate values.
  432.  
  433.                A numerical scheme for integrating the equations of motion
  434.                is presented, relevant coordinate  systems are defined and
  435.                algorithms for  transforming the various  input and output
  436.                parameters are summarized.
  437.  
  438.                All physical constants used by the software are included.
  439.  
  440.                The  algorithms  presented herein  support  the trajectory
  441.                mechanics  functions  in  the  program  source  code  with
  442.                examples that demonstrate their integrity.  A bibliography
  443.                is  provided  for  algorithm  extensions,  map  projection
  444.                related processes and routine mathematical processes.
  445.  
  446.  
  447.  
  448.  
  449.                Trajectory Maker                       Equations of Motion
  450.                Algorithm Supplement                                  2-1
  451.                ──────────────────────────────────────────────────────────
  452.                
  453.  
  454.  
  455.                2.  EQUATIONS OF MOTION
  456.  
  457.                The equations  of motion for  propagating trajectories are
  458.                modeled  in rectangular,  geocentric  inertial coordinates
  459.                with fundamental  plane coplanar with  the Earth's equator
  460.                and principle axis in the direction of the vernal equinox.
  461.  
  462.                Let x, y,  z denote these  coordinates with z-axis coinci-
  463.                dent with the Earth's  spin-axis, positive northwards; the
  464.                x-axis in the equatorial  plane, positive in the direction
  465.                of  the  vernal  equinox;  and  the  y-axis  lying  in the
  466.                equatorial  plane  in such  a  manner that  the  system is
  467.                right-handed.
  468.  
  469.                Further, let  i, j, k  be unit  vectors along the  x, y, z
  470.                axes, respectively.  Then, the position  of a point mass m
  471.                may be represented in vector form
  472.  
  473.                    r = xi + yj + zk
  474.  
  475.                and according  to  Newton's  second law,  the differential
  476.                equations of motion may be written
  477.  
  478.                      d2r
  479.                    m ─── =  Σ F
  480.                      dt2
  481.  
  482.                where Σ F is the sum  of the external forces acting on the
  483.                mass m.
  484.  
  485.                In this paper the forces  due to the Earth's gravitational
  486.                potential  function   and  its   atmosphere  are  treated.
  487.                Letting G and D denote these  forces per unit mass, we may
  488.                write
  489.  
  490.                    d2r
  491.                    ─── = G + D
  492.                    dt2
  493.  
  494.                By  substitution  of  variables,   this  equation  may  be
  495.                decomposed  into  the following  first  order differential
  496.                equations,
  497.  
  498.                    dr
  499.                    ── = v
  500.                    dt
  501.  
  502.  
  503.  
  504.  
  505.                Trajectory Maker                       Equations of Motion
  506.                Algorithm Supplement                                  2-2
  507.                ──────────────────────────────────────────────────────────
  508.                
  509.  
  510.  
  511.  
  512.                    dv
  513.                    ── = G + D
  514.                    dt
  515.  
  516.                where,  using  Newtonian   differentiation,  the  velocity
  517.                vector is
  518.                        .    .    .
  519.                    v = xi + yj + zk
  520.  
  521.                Expressed in component form, these two vector differential
  522.                equations become
  523.  
  524.                    dx   .            dy   .            dz   .
  525.                    ── = x ,          ── = y  ,         ── = z
  526.                    dt                dt                dt
  527.  
  528.                     .                 .                 .
  529.                    dx                dy                dz
  530.                    ── = Gx + Dx ,    ── = Gy + Dy ,    ── = Gz + Dz
  531.                    dt                dt                dt
  532.  
  533.                A  solution to  these  six differential  equations  may be
  534.                obtained,  once the  components  of the  gravity  and drag
  535.                acceleration components  are expressed  as known functions
  536.                of x, y,  z.  A time  series solution is  called an ephem-
  537.                eris.  The actual path of the mass m is its trajectory.
  538.  
  539.  
  540.  
  541.  
  542.                Trajectory Maker                             Gravity Model
  543.                Algorithm Supplement                                  3-1
  544.                ──────────────────────────────────────────────────────────
  545.                
  546.  
  547.  
  548.                3.  GRAVITY PERTURBATION MODEL
  549.  
  550.                We suppose that the Earth's  figure is an oblate spheroid,
  551.                and let P be point at a distance r from its center of mass
  552.                and let φ'  be the declination of  P.  Then, the potential
  553.                at P may  be expressed as a  series of spherical harmonics
  554.                of the form
  555.  
  556.                           µ ┌─     ∞                        ─┐
  557.                    V = - ── │ 1 -  Σ  Jn ( R/r)n  Pn(sin φ') │
  558.                           r └─    n=2                       ─┘
  559.  
  560.                where µ is the Earth's  gravitational mass constant, R its
  561.                equatorial radius, Jn the  zonal harmonic coefficients and
  562.                Pn the associated  Legendre polynomials of  order n.  This
  563.                form of  the Earth's  potential function  is known  as the
  564.                Vinti potential function [1].
  565.  
  566.                The gravitational  mass constant  and Earth  radii used in
  567.                this  paper are  from the  WGS-72  constants [3]  shown in
  568.                Table 3-1.
  569.  
  570.  
  571.                          Table 3-1.  WGS-72 Earth Constants
  572.  
  573.                  ╒═══╤═══════════════════════╤═════════════════════╕
  574.                  │   │      English          │      Metric         │
  575.                  ├───┼───────────────────────┼─────────────────────┤
  576.                  │ µ │ 62750.16789 nm3/sec2  │ 398600.50 km3/sec2  │
  577.                  ├───┼───────────────────────┼─────────────────────┤
  578.                  │ R │  3443.91739 nm        │   6378.135 km       │
  579.                  ├───┼───────────────────────┴─────────────────────┤
  580.                  │ Ω │  0.7292115147 (10)-4  rad/sec               │
  581.                  ├───┼─────────────────────────────────────────────┤
  582.                  │ f │  1/298.26  Dimensionless                    │
  583.                  └───┴─────────────────────────────────────────────┘
  584.  
  585.  
  586.                The  harmonic  coefficients are  the  Smithsonian  1973 II
  587.                coefficients [2] tabulated in Table 3-2.
  588.  
  589.  
  590.  
  591.  
  592.                Trajectory Maker                             Gravity Model
  593.                Algorithm Supplement                                  3-2
  594.                ──────────────────────────────────────────────────────────
  595.                
  596.  
  597.  
  598.                  Table 3-2.  Smithsonian Zonal Harmonic Coefficients
  599.  
  600.                  ╒════════════════════════╦════════════════════════╕
  601.                  │    J2 = 1082.636e-6    ║    J13 =  -0.336e-6    │
  602.                  │    J3 =   -2.540e-6    ║    J14 =   0.101e-6    │
  603.                  │    J4 =   -1.619e-6    ║    J15 =   0.104e-6    │
  604.                  │    J5 =   -0.230e-6    ║    J16 =   0.043e-6    │
  605.                  │    J6 =    0.552e-6    ║    J17 =  -0.227e-6    │
  606.                  │    J7 =   -0.345e-6    ║    J18 =  -0.077e-6    │
  607.                  │    J8 =   -0.204e-6    ║    J19 =   0.083e-6    │
  608.                  │    J9 =   -0.162e-6    ║    J20 =  -0.108e-6    │
  609.                  │    J10 =  -0.232e-6    ║    J21 =  -0.070e-6    │
  610.                  │    J11 =   0.317e-6    ║    J22 =   0.075e-6    │
  611.                  │    J12 =  -0.196e-6    ║    J23 =   0.111e-6    │
  612.                  └────────────────────────╨────────────────────────┘
  613.  
  614.  
  615.                Now, the distance  from the Earth's center  to the mass m,
  616.                coincident with the point P,  is the magnitude of r, given
  617.                by
  618.  
  619.                   r = √(x2 + y2 + z2)
  620.  
  621.                and the sine of its declination is
  622.  
  623.                    sin φ' = z/r
  624.  
  625.                Since  the gravitational  force  is conservative,  the net
  626.                gravitational  acceleration G  acting on  the mass  is the
  627.                negative gradient of V
  628.  
  629.                    G = -   V
  630.  
  631.                which can be decomposed, by differentiating, into the form
  632.  
  633.                    G = -(µ/r3)r -   U
  634.  
  635.                The first term in  this expression represents the dominant
  636.                attraction the  Earth would  exhibit if  it were  a sphere
  637.                comprised  of  homogeneous concentric  shells.   This term
  638.                accounts for the non-perturbed conic motion of m.
  639.  
  640.                The function U  is the disturbing  function whose gradient
  641.                represents  the  gravity  perturbation  acceleration which
  642.                causes  deviations  from  the  conic  motion  of  m.  This
  643.                function is given by
  644.  
  645.  
  646.  
  647.  
  648.                Trajectory Maker                             Gravity Model
  649.                Algorithm Supplement                                  3-3
  650.                ──────────────────────────────────────────────────────────
  651.                
  652.  
  653.  
  654.                           ∞
  655.                    U = µ  Σ  Jn Rn r-(n+1) Pn(sin φ')
  656.                          n=2
  657.  
  658.                It follows from the  gradient operator that the components
  659.                of the gravitational perturbation acceleration are
  660.  
  661.                    δGx ≡ -   U/  x = (µx)/r3 Σ Jn(R/r)n fn(φ')
  662.  
  663.                    δGy ≡ -   U/  y = (µy)/r3 Σ Jn(R/r)n fn(φ')
  664.  
  665.                    δGz ≡ -   U/  z = (µz)/r3 Σ Jn(R/r)n gn(φ')
  666.  
  667.                where
  668.  
  669.                 fn(sin φ') = (n + 1)Pn(sin φ') + sin φ Pn'(sin φ')
  670.  
  671.                 gn(sin φ') = (n + 1)Pn(sin φ') - cos φ cot φ Pn'(sin φ')
  672.  
  673.                Whittaker  and  Watson  [5]   derive  the  following  four
  674.                recurrence formulas for the complex variable u:
  675.  
  676.                    Pn+1'(u) - uPn'(u) = (n + 1)Pn(u)                  (1)
  677.  
  678.                    (u2 - 1)Pn'(u) = nuPn(u) - nPn-1(u)                (2)
  679.  
  680.                    (n + 1)Pn+1(u) - (2n + 1)uPn(u) + nPn-1(u) = 0     (3)
  681.  
  682.                    Pn+1'(u) - Pn-1'(u) = (2n + 1)Pn(u)                (4)
  683.  
  684.                Using the first recurrence formula yields
  685.  
  686.                    fn(sin φ') = Pn+1'(sin φ')
  687.  
  688.                and the second and third formulas yield
  689.  
  690.                    gn(sin φ') = (n + 1)Pn+1(sin φ')/sin φ'
  691.  
  692.                Therefore, the components of the perturbation acceleration
  693.                become
  694.  
  695.                    δGx = (µx)/r3 Σ Jn(R/r)n Pn+1'(sin φ')
  696.  
  697.                    δGy = (µy)/r3 Σ Jn(R/r)n Pn+1'(sin φ')
  698.  
  699.                    δGz = (µ)/r2  Σ Jn(R/r)n(n + 1)Pn+1(sin φ')
  700.  
  701.                In order to numerically  evaluate these expressions it is,
  702.  
  703.  
  704.  
  705.  
  706.                Trajectory Maker                             Gravity Model
  707.                Algorithm Supplement                                  3-4
  708.                ──────────────────────────────────────────────────────────
  709.                
  710.  
  711.  
  712.                of  course, necessary  to limit  the  series' to  a finite
  713.                number of  terms.  Suppose, then,  that the  degree of the
  714.                polynomials in R/r is specified, say  M ≥ 2.  Then, with a
  715.                shift of indices we have
  716.  
  717.                                         M-1
  718.                    δGx = (µx)/r3 (R/r)2  Σ  ai(R/r)i-1
  719.                                         i=1
  720.  
  721.                                         M-1
  722.                    δGy = (µy)/r3 (R/r)2  Σ  ai(R/r)i-1
  723.                                         i=1
  724.  
  725.                                         M-1
  726.                    δGz = (µ)/r2  (R/r)2  Σ  bi(R/r)i-1
  727.                                         i=1
  728.  
  729.                where
  730.  
  731.                    ai = Ji+1 Pi+2'(sin φ')
  732.  
  733.                    bi = Ji+1(i + 2)Pi+2(sin φ')
  734.  
  735.                The Legendre polynomials  and their derivatives  for n ≥ 3
  736.                may  be respectively  computed from  the third  and fourth
  737.                recurrence formulas once it is noted that [5]
  738.  
  739.                    P1(u) = u
  740.  
  741.                    P2(u) = 1/2 (3u2 - 1)
  742.  
  743.                and consequently
  744.  
  745.                    P1'(u) = 1
  746.  
  747.                    P2'(u) = 3u
  748.  
  749.                In order to illustrate the effects of various combinations
  750.                of gravity  harmonic coefficients,  we consider  a nominal
  751.                Defense Meteorological Satellite Program mission.
  752.  
  753.                Its state  vector, at  orbit insertion,  is represented in
  754.                the following Earth-centered-inertial coordinates:
  755.  
  756.  
  757.  
  758.  
  759.                Trajectory Maker                             Gravity Model
  760.                Algorithm Supplement                                  3-5
  761.                ──────────────────────────────────────────────────────────
  762.                
  763.  
  764.  
  765.                                             .
  766.                    x =   442.151588 nm      x =  0.512255 nm/sec
  767.                                             .
  768.                    y =  1387.396376 nm      y = -3.732104 nm/sec
  769.                                             .
  770.                    z = -3611.173591 nm      z = -1.373147 nm/sec
  771.  
  772.                The  satellite's  near circular  orbit  is sun-synchronous
  773.                with altitude of about 450 nm.
  774.  
  775.                The orbit ephemerides were propagated from orbit insertion
  776.                for a duration of 25000 seconds, in increments of 100 sec.
  777.                This step-size is optimal, in this case, since any smaller
  778.                step  does not  significantly  improve the  solution.  The
  779.                trajectory ephemeris spans eight orbits.
  780.  
  781.                Various gravity  models with reduced  orders (less coeffi-
  782.                cients)  were compared  with  the model  of  highest order
  783.                having the 22 coefficients listed in Table 2 above.
  784.  
  785.                The state vector  of the highest  order model was compared
  786.                with the lower order  models by differencing their respec-
  787.                tive  position   vector  components,   and  computing  the
  788.                magnitude  of  the resulting  difference  vector.   In the
  789.                figures the difference between a trajectory generated with
  790.                a gravity model of order M and one of order N is indicated
  791.                in the legend by the symbol DJMJN.
  792.  
  793.                The first group of the  figures, namely, Figures 3-1, 3-2,
  794.                3-3, 3-4, 3-5, 3-6, illustrate differences in ephemerides,
  795.                sequentially beginning with the  J2 coefficient and ending
  796.                with coefficients up to J23.   It is seen that this family
  797.                of curves has  a secular growth over  the time span, while
  798.                exhibiting a  tendency towards  convergence; however, that
  799.                is not to say that  any given order gravity model produces
  800.                more accurate trajectories than a lower order model, as
  801.                for  example, comparison  of the  sixth and  seventh order
  802.                models  in   Figure  3-2.   Figure   3-3  illustrates  the
  803.                beginning  of  the  loss of  data  integrity  as round-off
  804.                accumulates, becoming  more pronounced in  Figures 3-4 and
  805.                3-5.
  806.  
  807.                The scales in  the above group of  figures vary across the
  808.                figures by  two orders of  magnitude.  In  order to illus-
  809.                trate variations across scales, the next group of figures,
  810.                namely,  Figures 3-7,  3-8, 3-9,  and 3-10  were composed.
  811.                Not all differences are illustrated, but enough to get the
  812.                idea.
  813.  
  814.  
  815.  
  816.  
  817.                Trajectory Maker                             Gravity Model
  818.                Algorithm Supplement                                  3-6
  819.                ──────────────────────────────────────────────────────────
  820.                
  821.  
  822.  
  823.                The  gravitational effects  on  other orbit  parameters is
  824.                given in detail in Reynolds [4].
  825.  
  826.  
  827.  
  828.  
  829.                Trajectory Maker                             Gravity Model
  830.                Algorithm Supplement                                  3-7
  831.                ──────────────────────────────────────────────────────────
  832.                
  833.  
  834.  
  835.  
  836.  
  837.  
  838.  
  839.  
  840.  
  841.  
  842.  
  843.  
  844.  
  845.  
  846.  
  847.  
  848.  
  849.  
  850.  
  851.  
  852.  
  853.  
  854.  
  855.  
  856.  
  857.  
  858.  
  859.  
  860.  
  861.  
  862.  
  863.  
  864.  
  865.  
  866.  
  867.  
  868.  
  869.  
  870.  
  871.  
  872.  
  873.  
  874.  
  875.  
  876.  
  877.  
  878.  
  879.                Figure 3-1.  Magnitude of Vector Differences Between the
  880.                             23rd and 2nd, 3rd, 4th and 5th Order Gravity
  881.                             Models.
  882.  
  883.  
  884.  
  885.  
  886.                Trajectory Maker                             Gravity Model
  887.                Algorithm Supplement                                  3-8
  888.                ──────────────────────────────────────────────────────────
  889.                
  890.  
  891.  
  892.  
  893.  
  894.  
  895.  
  896.  
  897.  
  898.  
  899.  
  900.  
  901.  
  902.  
  903.  
  904.  
  905.  
  906.  
  907.  
  908.  
  909.  
  910.  
  911.  
  912.  
  913.  
  914.  
  915.  
  916.  
  917.  
  918.  
  919.  
  920.  
  921.  
  922.  
  923.  
  924.  
  925.  
  926.  
  927.  
  928.  
  929.  
  930.  
  931.  
  932.  
  933.  
  934.  
  935.  
  936.                Figure 3-2.  Magnitude of Vector Differences Between the
  937.                             23rd and 6th, 7th, 8th and 9th Order Gravity
  938.                             Models.
  939.  
  940.  
  941.  
  942.  
  943.                Trajectory Maker                             Gravity Model
  944.                Algorithm Supplement                                  3-9
  945.                ──────────────────────────────────────────────────────────
  946.                
  947.  
  948.  
  949.  
  950.  
  951.  
  952.  
  953.  
  954.  
  955.  
  956.  
  957.  
  958.  
  959.  
  960.  
  961.  
  962.  
  963.  
  964.  
  965.  
  966.  
  967.  
  968.  
  969.  
  970.  
  971.  
  972.  
  973.  
  974.  
  975.  
  976.  
  977.  
  978.  
  979.  
  980.  
  981.  
  982.  
  983.  
  984.  
  985.  
  986.  
  987.  
  988.  
  989.  
  990.  
  991.  
  992.  
  993.                Figure 3-3.  Magnitude of Vector Differences Between the
  994.                             23rd and 10th, 11th, 12th and 13th Order
  995.                             Gravity Models.
  996.  
  997.  
  998.  
  999.  
  1000.                Trajectory Maker                             Gravity Model
  1001.                Algorithm Supplement                                  3-10
  1002.                ──────────────────────────────────────────────────────────
  1003.                
  1004.  
  1005.  
  1006.  
  1007.  
  1008.  
  1009.  
  1010.  
  1011.  
  1012.  
  1013.  
  1014.  
  1015.  
  1016.  
  1017.  
  1018.  
  1019.  
  1020.  
  1021.  
  1022.  
  1023.  
  1024.  
  1025.  
  1026.  
  1027.  
  1028.  
  1029.  
  1030.  
  1031.  
  1032.  
  1033.  
  1034.  
  1035.  
  1036.  
  1037.  
  1038.  
  1039.  
  1040.  
  1041.  
  1042.  
  1043.                
  1044.  
  1045.  
  1046.  
  1047.  
  1048.  
  1049.  
  1050.                Figure 3-4.   Magnitude of Vector  Differences Between the
  1051.                             23rd  and  14th, 15th,  16th  and  17th Order
  1052.                             Gravity Models.
  1053.  
  1054.  
  1055.  
  1056.  
  1057.                Trajectory Maker                             Gravity Model
  1058.                Algorithm Supplement                                  3-11
  1059.                ──────────────────────────────────────────────────────────
  1060.                
  1061.  
  1062.  
  1063.  
  1064.  
  1065.  
  1066.  
  1067.  
  1068.  
  1069.  
  1070.  
  1071.  
  1072.  
  1073.  
  1074.  
  1075.  
  1076.  
  1077.  
  1078.  
  1079.  
  1080.  
  1081.  
  1082.  
  1083.  
  1084.  
  1085.  
  1086.  
  1087.  
  1088.  
  1089.  
  1090.  
  1091.  
  1092.  
  1093.  
  1094.  
  1095.  
  1096.  
  1097.  
  1098.  
  1099.                
  1100.                
  1101.  
  1102.  
  1103.  
  1104.  
  1105.  
  1106.                
  1107.                Figure 3-5.   Magnitude of Vector  Differences Between the
  1108.                             23rd  and  18th, 19th,  20th  and  21st Order
  1109.                             Gravity Models.
  1110.  
  1111.  
  1112.  
  1113.  
  1114.                Trajectory Maker                             Gravity Model
  1115.                Algorithm Supplement                                  3-12
  1116.                ──────────────────────────────────────────────────────────
  1117.                
  1118.  
  1119.  
  1120.  
  1121.  
  1122.  
  1123.  
  1124.  
  1125.  
  1126.  
  1127.  
  1128.  
  1129.  
  1130.  
  1131.  
  1132.  
  1133.  
  1134.  
  1135.  
  1136.  
  1137.  
  1138.  
  1139.  
  1140.  
  1141.  
  1142.  
  1143.  
  1144.  
  1145.  
  1146.  
  1147.  
  1148.  
  1149.  
  1150.  
  1151.  
  1152.                
  1153.                
  1154.  
  1155.  
  1156.  
  1157.  
  1158.  
  1159.  
  1160.  
  1161.                
  1162.                
  1163.                
  1164.                Figure 3-6.   Magnitude of Vector  Differences Between the
  1165.                             23rd  and  19th, 20th,  21st  and  22nd Order
  1166.                             Gravity Models.
  1167.  
  1168.  
  1169.  
  1170.  
  1171.                Trajectory Maker                             Gravity Model
  1172.                Algorithm Supplement                                  3-13
  1173.                ──────────────────────────────────────────────────────────
  1174.                
  1175.  
  1176.  
  1177.  
  1178.  
  1179.  
  1180.  
  1181.  
  1182.  
  1183.  
  1184.  
  1185.  
  1186.  
  1187.  
  1188.  
  1189.  
  1190.  
  1191.  
  1192.  
  1193.  
  1194.  
  1195.  
  1196.  
  1197.  
  1198.  
  1199.  
  1200.  
  1201.  
  1202.  
  1203.  
  1204.  
  1205.  
  1206.  
  1207.  
  1208.  
  1209.  
  1210.  
  1211.  
  1212.  
  1213.  
  1214.  
  1215.  
  1216.  
  1217.  
  1218.                
  1219.                
  1220.                
  1221.                Figure 3-7.   Magnitude of Vector  Differences Between the
  1222.                             23rd and 2nd,  4th, 6th and  8th Order (Even)
  1223.                             Gravity Models.
  1224.  
  1225.  
  1226.  
  1227.  
  1228.                Trajectory Maker                             Gravity Model
  1229.                Algorithm Supplement                                  3-14
  1230.                ──────────────────────────────────────────────────────────
  1231.                
  1232.  
  1233.  
  1234.  
  1235.  
  1236.  
  1237.  
  1238.  
  1239.  
  1240.  
  1241.  
  1242.  
  1243.  
  1244.  
  1245.  
  1246.  
  1247.  
  1248.  
  1249.  
  1250.  
  1251.  
  1252.  
  1253.  
  1254.  
  1255.  
  1256.  
  1257.  
  1258.  
  1259.  
  1260.  
  1261.  
  1262.  
  1263.  
  1264.  
  1265.  
  1266.  
  1267.  
  1268.  
  1269.  
  1270.  
  1271.  
  1272.  
  1273.  
  1274.  
  1275.  
  1276.  
  1277.  
  1278.                Figure 3-8.   Magnitude of Vector  Differences Between the
  1279.                             23rd and  3rd, 5th,  7th and  9th Order (Odd)
  1280.                             Gravity Models.
  1281.  
  1282.  
  1283.  
  1284.  
  1285.                Trajectory Maker                             Gravity Model
  1286.                Algorithm Supplement                                  3-15
  1287.                ──────────────────────────────────────────────────────────
  1288.                
  1289.  
  1290.  
  1291.  
  1292.  
  1293.  
  1294.  
  1295.  
  1296.  
  1297.  
  1298.  
  1299.  
  1300.  
  1301.  
  1302.  
  1303.  
  1304.  
  1305.  
  1306.  
  1307.  
  1308.  
  1309.  
  1310.  
  1311.  
  1312.  
  1313.  
  1314.  
  1315.  
  1316.  
  1317.  
  1318.  
  1319.  
  1320.  
  1321.  
  1322.  
  1323.  
  1324.  
  1325.  
  1326.  
  1327.  
  1328.  
  1329.  
  1330.  
  1331.  
  1332.  
  1333.  
  1334.  
  1335.                Figure 3-9.   Magnitude of Vector  Differences Between the
  1336.                             23rd  and  2nd,  6th,  10th  and  14th  Order
  1337.                             Gravity Models.
  1338.  
  1339.  
  1340.  
  1341.                Trajectory Maker                             Gravity Model
  1342.                Algorithm Supplement                                  3-16
  1343.                ──────────────────────────────────────────────────────────
  1344.                
  1345.  
  1346.  
  1347.  
  1348.  
  1349.  
  1350.  
  1351.  
  1352.  
  1353.  
  1354.  
  1355.  
  1356.  
  1357.  
  1358.  
  1359.  
  1360.  
  1361.  
  1362.  
  1363.  
  1364.  
  1365.  
  1366.  
  1367.  
  1368.  
  1369.  
  1370.  
  1371.  
  1372.  
  1373.  
  1374.  
  1375.  
  1376.  
  1377.  
  1378.  
  1379.  
  1380.  
  1381.  
  1382.  
  1383.  
  1384.  
  1385.  
  1386.  
  1387.  
  1388.  
  1389.  
  1390.  
  1391.                Figure 3-10.  Magnitude of  Vector Differences Between the
  1392.                              23rd  and  5th, 10th,  15th  and  20th Order
  1393.                              Gravity Models.
  1394.  
  1395.  
  1396.  
  1397.  
  1398.                Trajectory Maker                             Gravity Model
  1399.                Algorithm Supplement                                  3-17
  1400.                ──────────────────────────────────────────────────────────
  1401.                
  1402.  
  1403.  
  1404.                                        References
  1405.  
  1406.                1.  Anon., Space Flight Handbooks, Volume I, NASA SP-33
  1407.                    Part 1, Office of Scientific and Technical Inform-
  1408.                    ation, NASA, Washington D. C., 1963.
  1409.  
  1410.                2.  Gaposchkin, F. M., ed., 1973 Smithsonian Standard
  1411.                    Earth (III), Smithsonian Astrophysical Observatory,
  1412.                    Special Report 353, Cambridge, Ma., November 28, 1973.
  1413.  
  1414.                3.  Hieb, H., Western Test Range Geodetic Coordinate
  1415.                    Manual, Part 1, Federal Electric Corporation, Feb.
  1416.                    1980.
  1417.  
  1418.                4.  Reynolds, H. B., High-Order Gravitational Perturbation
  1419.                    Effects On DMSP Extended Mission Ephemerides, OAO
  1420.                    Corporation, Technical Operating Report, El Segundo,
  1421.                    Ca., Jan. 3, 1985.
  1422.  
  1423.                5.  Whittaker, F. and G. Watson, A Course in Modern
  1424.                    Analysis, Cambridge University Press, Fourth Ed. 1963.
  1425.                
  1426.  
  1427.  
  1428.  
  1429.                Trajectory Maker                                Drag Model
  1430.                Algorithm Supplement                                  4-1
  1431.                ──────────────────────────────────────────────────────────
  1432.                
  1433.  
  1434.  
  1435.                4.  DRAG PERTURBATION MODEL
  1436.  
  1437.                The  aerodynamic drag  force  acting on  an  object varies
  1438.                jointly with  the atmospheric  density σ,  the object drag
  1439.                coefficient CD,  the projected  cross sectional  area A of
  1440.                the object (normal to the  velocity vector), the square of
  1441.                its  speed relative  to  the rotating  atmosphere  va, and
  1442.                inversely  with its  mass  m.  The  drag  force acts  in a
  1443.                direction opposite that of the relative velocity vector.
  1444.  
  1445.                In vector form this force per unit mass may be written
  1446.  
  1447.                            CD A
  1448.                    D = - σ ────  va va
  1449.                             2m
  1450.  
  1451.                justifiably, since the object  imparts a velocity of order
  1452.                va  to a  mass of  air Aσva  in each  second; the  rate of
  1453.                change of momentum, equal to the drag force on the object,
  1454.                is thus Aσva².  The empirical drag coefficient is based on
  1455.                such  things  as  the object's  shape,  altitude  and Mach
  1456.                number.
  1457.  
  1458.                Taking into account local  wind conditions with the vector
  1459.                q, the velocity relative to the rotating air mass is
  1460.  
  1461.                    va = v - Ω x r - q
  1462.  
  1463.                where Ω  = Ωk is  the Earth's rotation  rate vector having
  1464.                magnitude 0.7292115147 (10)-4  rad/sec.
  1465.  
  1466.                It is usual to define a ballistic coefficient ß by
  1467.  
  1468.                         m
  1469.                    ß = ────
  1470.                        CD A
  1471.  
  1472.                and write the drag force as
  1473.  
  1474.                            1
  1475.                    D = -  ─── σ va va
  1476.                           2 ß
  1477.  
  1478.                which in component form,  after expanding the above vector
  1479.                product, becomes
  1480.  
  1481.  
  1482.  
  1483.  
  1484.                Trajectory Maker                                Drag Model
  1485.                Algorithm Supplement                                  4-2
  1486.                ──────────────────────────────────────────────────────────
  1487.                
  1488.  
  1489.  
  1490.                               .
  1491.                    Dx = -σva[(x + Ωy) - q1]/(2ß)
  1492.                               .
  1493.                    Dy = -σva[(y - Ωx) - q2]/(2ß)
  1494.                               .
  1495.                    Dz = -σva[ z - q3 ]/(2ß)
  1496.  
  1497.                where
  1498.                            .                  .                 .
  1499.                   va = √{[(x + Ωy) - q1]2 + [(y - Ωx) - q2]2 + (z - q3)2}
  1500.  
  1501.                and  q1,  q2,  q3  denote  the  wind  velocity  components
  1502.                relative to the rotating atmosphere.
  1503.  
  1504.                Ballistic  coefficients  have  an  intuitive  appeal since
  1505.                their units,  measured in  kg/m²,  resemble density units.
  1506.                When working with English units, the ballistic coefficient
  1507.                is usually expressed in mean  sea level weight rather than
  1508.                mass, that is lb/ft².   The following are typical examples
  1509.                that may be encountered:
  1510.  
  1511.                    o A space based radar antenna: 0.144, to 0.448 lb/ft²
  1512.                    o Space Shuttle Discovery: 42 lb/ft²
  1513.                    o Ballistic missiles: 25, 75, 2150 lb/ft²
  1514.  
  1515.                When  computing  a  ballistic  coefficients  for  orbiting
  1516.                objects,  the  parameter  with  the  most  uncertainty  is
  1517.                usually the drag coefficient.  Although this is not always
  1518.                the  case.  For  example, Reynolds  [2] considers  a space
  1519.                based radar system  having a rotating  antenna whose shape
  1520.                is an arbitrary surface of revolution.
  1521.  
  1522.                Wertz  [5] suggests  that  in the  absence  of aposteriori
  1523.                knowledge of the  drag coefficient that a  value of 2.0 be
  1524.                used.   Waison [4]  reports that  it should  be consistent
  1525.                with the value  used to derive  the atmospheric model, and
  1526.                that for the 1959  ARDC model one should  use 2.0; for the
  1527.                1964 Jacchia model, a value of 2.2 should be used.
  1528.  
  1529.                For reentry vehicles, the drag coefficient may be provided
  1530.                as  a function  of Mach  number  M, as  in the  Table 4-1,
  1531.                together with the vehicle's mass and effective area.
  1532.  
  1533.  
  1534.  
  1535.  
  1536.                Trajectory Maker                                Drag Model
  1537.                Algorithm Supplement                                  4-3
  1538.                ──────────────────────────────────────────────────────────
  1539.                
  1540.  
  1541.  
  1542.                      Table 4-1.  Drag Coefficient vs. Mach Number
  1543.  
  1544.                         Mass = 35.8365 slugs,   Area = 19.8 ft²
  1545.                 ┌────╥─────┬─────┬─────┬─────┬─────┬──────┬─────┬─────┐
  1546.                 │ M  ║ 0   │ .25 │ .50 │ .75 │ 1.0 │ 1.25 │ 1.5 │ 100 │
  1547.                 ├────╫─────┼─────┼─────┼─────┼─────┼──────┼─────┼─────┤
  1548.                 │ CD ║ .38 │ .40 │ .44 │ .55 │ .72 │ .76  │ .77 │ .77 │
  1549.                 └────╨─────┴─────┴─────┴─────┴─────┴──────┴─────┴─────┘
  1550.  
  1551.  
  1552.                Mach number  is function of  the speed of  sound, which in
  1553.                turn is a function of temperature.  The expression adopted
  1554.                for the speed of sound in Reference [1] is
  1555.  
  1556.                    Cs = √(  R* TM/M0)
  1557.  
  1558.                where   is the  ratio of specific heat  of air at constant
  1559.                pressure to that  at constant volume, and  taken to be 1.4
  1560.                (dimensionless) exactly, R*  the gas constant  taken to be
  1561.                8314.32  N.m/(kmol.K),  M0  the  sea-level  mean molecular
  1562.                weight  of  the  mixture of  gasses  taken  to  be 28.9644
  1563.                kg/kmol,  and  TM  is  the  molecular  scale  temperature,
  1564.                modeled as a linear, segmental arc.
  1565.  
  1566.                Putting all these things together results in
  1567.  
  1568.                    Cs = 20.04680276 √TM    (meters/sec)
  1569.  
  1570.                and the Mach number is
  1571.  
  1572.                    M = va/Cs
  1573.  
  1574.                The  concept of  speed of  sound progressively  looses its
  1575.                range  of applicability  at  high altitudes  and  for this
  1576.                reason the values  of speed of sound  are not computed for
  1577.                heights above 86 km.
  1578.  
  1579.                Table 4-2 is extracted from Table  9 of the U. S. Standard
  1580.                1976 abbreviated tables.
  1581.  
  1582.  
  1583.  
  1584.  
  1585.                Trajectory Maker                                Drag Model
  1586.                Algorithm Supplement                                  4-4
  1587.                ──────────────────────────────────────────────────────────
  1588.                
  1589.  
  1590.  
  1591.                                   Table 4-2
  1592.  
  1593.                            Molecular Scale Temperature
  1594.  
  1595.                       ╒═══════════════╤══════════════════╕
  1596.                       │ Geometric     │ Molecular Scale  │
  1597.                       │ Altitude (km) │ Temperature (K)  │
  1598.                       ├───────────────┼──────────────────┤
  1599.                       │     0.0000    │     288.150      │
  1600.                       │    11.0190    │     216.650      │
  1601.                       │    20.0631    │     216.650      │
  1602.                       │    32.1619    │     288.650      │
  1603.                       │    47.3500    │     270.650      │
  1604.                       │    51.4124    │     270.650      │
  1605.                       │    71.8019    │     214.650      │
  1606.                       │    86.0000    │     186.946      │
  1607.                       └───────────────┴──────────────────┘
  1608.  
  1609.  
  1610.                Since the U. S. STandard  models the temperature as linear
  1611.                segments, it is only  necessary to interpolate linearly in
  1612.                this  table  to  obtain  TM  for  arbitrary  altitude, and
  1613.                substitute the result  into the above  equation to compute
  1614.                the speed of sound.
  1615.  
  1616.  
  1617.                4.1  Atmospheric Density Profile
  1618.  
  1619.                Another complication  is that  density varies  in a fairly
  1620.                irregular  manner  with factors  such  as  solar activity,
  1621.                yearly  seasons, time  of  day, geographical  location and
  1622.                altitude.   Fortunately,   extensive  research   has  been
  1623.                performed  by various  agencies to  define models  for the
  1624.                principle atmospheric parameters.
  1625.  
  1626.                The 1976  U. S. Standard  Atmosphere [1]  is an idealized,
  1627.                steady state representation of the Earth's atmosphere from
  1628.                its surface  to 1000 km,  as it  is assumed to  exist in a
  1629.                period of moderate solar activity.  It consists of compre-
  1630.                hensive  tables  (approximately  1000  points)  of various
  1631.                atmospheric parameters, including  density, tabulated as a
  1632.                function of altitude.  The  density profile is illustrated
  1633.                in Figure 4-1 on a common logarithm scale.
  1634.  
  1635.                Since 1000 pairs is a lot  of data to handle, subsets were
  1636.                extracted and modeled  with the goal  of obtaining a rela-
  1637.                tively small subset that would yield sufficient accuracy.
  1638.  
  1639.  
  1640.  
  1641.  
  1642.                Trajectory Maker                                Drag Model
  1643.                Algorithm Supplement                                  4-5
  1644.                ──────────────────────────────────────────────────────────
  1645.                
  1646.  
  1647.  
  1648.                In effort to obtain  the smoothest possible representation
  1649.                for density, cubic spline functions were used to model the
  1650.                common logarithm of the density profile, with the realiza-
  1651.                tion  that,   during  orbit   propagation,  density  would
  1652.                necessarily  be  obtained from  the  logarithm  profile by
  1653.                exponentiation.
  1654.  
  1655.                Several  subsets were  extracted  from the  Standard Atmo-
  1656.                sphere, starting  with the  abbreviated tables, consisting
  1657.                of  21 pairs,  from the   Standard Atmosphere  [1].  These
  1658.                points are illustrated in Figure 4-2.
  1659.  
  1660.                The final set selected for the model, herein, was obtained
  1661.                by judiciously  selecting points from  the idealized model
  1662.                until the percent relative error in density was reduced to
  1663.                less than 1%.  The resulting set consisted of 48 altitude-
  1664.                density pairs illustrated in  Figure 4-3.  Note that these
  1665.                data are not equally spaced.  However,  besides being very
  1666.                fast, Thompson's  algorithm [3]  accounts for  this situa-
  1667.                tion.
  1668.  
  1669.                Figure 4-4 illustrates the percent relative error obtained
  1670.                when  using  these  two  density  profiles  to  model  the
  1671.                idealized atmosphere.  The dashed  curve shows that in the
  1672.                critical altitude  region (less  than 200  km) for reentry
  1673.                trajectories, the 21 point set  produces errors as high as
  1674.                6%.  On the  other hand, for  the 48 point  set the errors
  1675.                are less  than 1%.  This  is shown more  clearly in Figure
  1676.                4-5.
  1677.  
  1678.                Figure  4-6 compares  linear interpolation  (dashed curve)
  1679.                with spline interpolation (solid  curve).  It is seen that
  1680.                linear interpolation  produces errors  as high  as 4%.  It
  1681.                should be emphasized  that these error  curves are density
  1682.                errors rather than the logarithm of the density errors.
  1683.  
  1684.                Figure 4-7  isolates the error  for the 48-point  set on a
  1685.                magnified scale.   The figure shows  that the  error to be
  1686.                less than 1% on the entire altitude range.  It is this set
  1687.                that  comprises  the  density  profile  used  herein.  The
  1688.                actual values for this set are tabulated in Table 4-3.
  1689.  
  1690.                This study also found that the endpoint derivatives played
  1691.                a  significant  role  towards  achieving  the  1% accuracy
  1692.                criterion.   Thus, rather  than computing  the derivatives
  1693.                with  values from  Table  4-3, the  following  values were
  1694.                computed  using  a  finer  resolution  from  data  in  the
  1695.                comprehensive table:
  1696.  
  1697.  
  1698.  
  1699.  
  1700.                Trajectory Maker                                Drag Model
  1701.                Algorithm Supplement                                  4-6
  1702.                ──────────────────────────────────────────────────────────
  1703.                
  1704.  
  1705.  
  1706.  
  1707.                    (log σ(0))' = -0.041934
  1708.  
  1709.                    (log σ(1000))' = -0.001834
  1710.  
  1711.                where we  recall that the  splines were fit  to the common
  1712.                logarithm of the density values.
  1713.  
  1714.  
  1715.                        Table 4-3.  Atmospheric Density Profile
  1716.  
  1717.                 ╒════════════════════════╦═════════════════════════╕
  1718.                 │ Altitude   Density     ║ Altitude   Density      │
  1719.                 │  (km)      (kg/m3)     ║ (km)       (kg/m3)      │
  1720.                 ├────────────────────────╫─────────────────────────┤
  1721.                 │   0       1.2250       ║    90      3.416e-6     │
  1722.                 │   2       1.0066       ║   100      5.604e-7     │
  1723.                 │   4       8.1935e-1    ║   110      9.708e-8     │
  1724.                 │   6       6.6011e-1    ║   120      2.222e-8     │
  1725.                 │   8       5.2579e-1    ║   130      8.152e-9     │
  1726.                 │   10      4.1351e-1    ║   140      3.831e-9     │
  1727.                 │   12      3.1194e-1    ║   150      2.076e-9     │
  1728.                 │   14      2.2786e-1    ║   160      1.233e-9     │
  1729.                 │   16      1.6647e-1    ║   170      7.815e-10    │
  1730.                 │   18      1.2165e-1    ║   180      5.194e-10    │
  1731.                 │   20      8.8910e-2    ║   190      3.581e-10    │
  1732.                 │   25      4.0084e-2    ║   200      2.541e-10    │
  1733.                 │   30      1.8410e-2    ║   220      1.367e-10    │
  1734.                 │   35      8.4634e-3    ║   240      7.858e-11    │
  1735.                 │   40      3.9957e-3    ║   260      4.742e-11    │
  1736.                 │   45      1.9663e-3    ║   280      2.971e-11    │
  1737.                 │   50      1.0269e-3    ║   300      1.916e-11    │
  1738.                 │   45      1.9663e-3    ║   400      2.802e-12    │
  1739.                 │   60      3.0968e-4    ║   500      5.215e-13    │
  1740.                 │   65      1.6321e-4    ║   600      1.137e-13    │
  1741.                 │   70      8.2829e-5    ║   700      3.069e-14    │
  1742.                 │   75      3.9921e-5    ║   800      1.136e-14    │
  1743.                 │   80      1.8458e-5    ║   900      5.759e-15    │
  1744.                 │   85      8.2196e-6    ║  1000      3.561e-15    │
  1745.                 └────────────────────────╨─────────────────────────┘
  1746.  
  1747.  
  1748.  
  1749.  
  1750.                Trajectory Maker                                Drag Model
  1751.                Algorithm Supplement                                  4-7
  1752.                ──────────────────────────────────────────────────────────
  1753.                
  1754.  
  1755.  
  1756.  
  1757.  
  1758.  
  1759.  
  1760.  
  1761.  
  1762.  
  1763.  
  1764.  
  1765.  
  1766.  
  1767.  
  1768.  
  1769.  
  1770.  
  1771.  
  1772.  
  1773.  
  1774.  
  1775.  
  1776.  
  1777.  
  1778.  
  1779.  
  1780.  
  1781.  
  1782.  
  1783.  
  1784.  
  1785.  
  1786.  
  1787.  
  1788.  
  1789.  
  1790.  
  1791.  
  1792.  
  1793.  
  1794.  
  1795.  
  1796.  
  1797.  
  1798.  
  1799.  
  1800.  
  1801.  
  1802.                Fig. 4-1.  1976 U. S. Standard Atmosphere Density Profile.
  1803.  
  1804.  
  1805.  
  1806.  
  1807.                Trajectory Maker                                Drag Model
  1808.                Algorithm Supplement                                  4-8
  1809.                ──────────────────────────────────────────────────────────
  1810.                
  1811.  
  1812.  
  1813.  
  1814.  
  1815.  
  1816.  
  1817.  
  1818.  
  1819.  
  1820.  
  1821.  
  1822.  
  1823.  
  1824.  
  1825.  
  1826.  
  1827.  
  1828.  
  1829.  
  1830.  
  1831.  
  1832.  
  1833.  
  1834.  
  1835.  
  1836.  
  1837.  
  1838.  
  1839.  
  1840.  
  1841.  
  1842.  
  1843.  
  1844.  
  1845.  
  1846.  
  1847.  
  1848.  
  1849.  
  1850.  
  1851.  
  1852.  
  1853.  
  1854.                
  1855.  
  1856.                
  1857.  
  1858.                Fig. 4-2.  Abbreviated 1976 U. S. Standard Atmosphere
  1859.                           Density Profile.
  1860.  
  1861.  
  1862.  
  1863.  
  1864.                Trajectory Maker                                Drag Model
  1865.                Algorithm Supplement                                  4-9
  1866.                ──────────────────────────────────────────────────────────
  1867.                
  1868.  
  1869.  
  1870.  
  1871.  
  1872.  
  1873.  
  1874.  
  1875.  
  1876.  
  1877.  
  1878.  
  1879.  
  1880.  
  1881.  
  1882.  
  1883.  
  1884.  
  1885.  
  1886.  
  1887.  
  1888.  
  1889.  
  1890.  
  1891.  
  1892.  
  1893.  
  1894.  
  1895.  
  1896.  
  1897.  
  1898.  
  1899.  
  1900.  
  1901.  
  1902.  
  1903.  
  1904.  
  1905.  
  1906.  
  1907.  
  1908.  
  1909.  
  1910.  
  1911.  
  1912.  
  1913.                
  1914.                
  1915.                
  1916.                Fig. 4-3.  Trajectory Maker Density Profile Model.
  1917.  
  1918.  
  1919.  
  1920.  
  1921.                Trajectory Maker                                Drag Model
  1922.                Algorithm Supplement                                  4-10
  1923.                ──────────────────────────────────────────────────────────
  1924.                
  1925.  
  1926.  
  1927.  
  1928.  
  1929.  
  1930.  
  1931.  
  1932.  
  1933.  
  1934.  
  1935.  
  1936.  
  1937.  
  1938.  
  1939.  
  1940.  
  1941.  
  1942.  
  1943.  
  1944.  
  1945.  
  1946.  
  1947.  
  1948.  
  1949.  
  1950.  
  1951.  
  1952.  
  1953.  
  1954.  
  1955.  
  1956.  
  1957.  
  1958.  
  1959.  
  1960.  
  1961.  
  1962.  
  1963.  
  1964.  
  1965.  
  1966.  
  1967.  
  1968.  
  1969.  
  1970.  
  1971.  
  1972.                Fig. 4-4.  Cubic Spline Interpolation Error for the 21 and
  1973.                           48 Point Tables.
  1974.  
  1975.  
  1976.  
  1977.  
  1978.                Trajectory Maker                                Drag Model
  1979.                Algorithm Supplement                                  4-11
  1980.                ──────────────────────────────────────────────────────────
  1981.                
  1982.  
  1983.  
  1984.  
  1985.  
  1986.  
  1987.  
  1988.  
  1989.  
  1990.  
  1991.  
  1992.  
  1993.  
  1994.  
  1995.  
  1996.  
  1997.  
  1998.  
  1999.  
  2000.  
  2001.  
  2002.  
  2003.  
  2004.  
  2005.  
  2006.  
  2007.  
  2008.  
  2009.  
  2010.  
  2011.  
  2012.  
  2013.  
  2014.  
  2015.  
  2016.  
  2017.  
  2018.  
  2019.  
  2020.  
  2021.  
  2022.  
  2023.  
  2024.  
  2025.                
  2026.  
  2027.                
  2028.  
  2029.                Fig. 4-5.  Cubic Spline Interpolation Error for the 21 and
  2030.                           48 Point Tables.
  2031.  
  2032.  
  2033.  
  2034.  
  2035.                Trajectory Maker                                Drag Model
  2036.                Algorithm Supplement                                  4-12
  2037.                ──────────────────────────────────────────────────────────
  2038.                
  2039.  
  2040.  
  2041.  
  2042.  
  2043.  
  2044.  
  2045.  
  2046.  
  2047.  
  2048.  
  2049.  
  2050.  
  2051.  
  2052.  
  2053.  
  2054.  
  2055.  
  2056.  
  2057.  
  2058.  
  2059.  
  2060.  
  2061.  
  2062.  
  2063.  
  2064.  
  2065.  
  2066.  
  2067.  
  2068.  
  2069.  
  2070.  
  2071.  
  2072.  
  2073.  
  2074.  
  2075.  
  2076.  
  2077.  
  2078.  
  2079.  
  2080.  
  2081.  
  2082.  
  2083.                
  2084.  
  2085.  
  2086.                Fig. 4-6.  Cubic Spline and Linear Interpolation Error
  2087.                           for the 48 Point Table.
  2088.  
  2089.  
  2090.  
  2091.  
  2092.                Trajectory Maker                                Drag Model
  2093.                Algorithm Supplement                                  4-13
  2094.                ──────────────────────────────────────────────────────────
  2095.                
  2096.  
  2097.  
  2098.  
  2099.  
  2100.  
  2101.  
  2102.  
  2103.  
  2104.  
  2105.  
  2106.  
  2107.  
  2108.  
  2109.  
  2110.  
  2111.  
  2112.  
  2113.  
  2114.  
  2115.  
  2116.  
  2117.  
  2118.  
  2119.  
  2120.  
  2121.  
  2122.  
  2123.  
  2124.  
  2125.  
  2126.  
  2127.  
  2128.  
  2129.  
  2130.  
  2131.  
  2132.  
  2133.  
  2134.  
  2135.  
  2136.  
  2137.  
  2138.                
  2139.  
  2140.                
  2141.  
  2142.                
  2143.                Fig.  4-7.  Cubic  Spline Interpolation  Error for  the 48
  2144.                           Point Table.
  2145.  
  2146.  
  2147.  
  2148.  
  2149.                Trajectory Maker                                Drag Model
  2150.                Algorithm Supplement                                  4-14
  2151.                ──────────────────────────────────────────────────────────
  2152.                
  2153.  
  2154.  
  2155.                4.2  Wind Profile
  2156.  
  2157.                Local atmospheric winds are typically supplied in the form
  2158.                of a table which contains the wind velocity magnitude q(h)
  2159.                and direction Θ(h)  as functions of  altitude.  The direc-
  2160.                tion of the wind is defined as the true bearing from which
  2161.                it is blowing.   The east, north and  up components of the
  2162.                wind velocity are, respectively,
  2163.  
  2164.                    qe = -q(h) sin Θ(h)
  2165.  
  2166.                    qn = -q(h) cos Θ(h)
  2167.  
  2168.                    qu = 0
  2169.  
  2170.                For the purpose of illustration,  a couple of typical wind
  2171.                profiles from the Lompoc region are illustrated in Figures
  2172.                4-8  and   4-9.   The  symbols   represent  observed  wind
  2173.                velocities that were mapped into east and north components
  2174.                by the above equations.   Cubic spline functions were used
  2175.                to model these  data and several  values were computed and
  2176.                plotted at  intermediate values  of altitude.   Two things
  2177.                are noteworthy - that  splines provide a smooth represent-
  2178.                ation of the data, and that Lompoc is a windy place.
  2179.  
  2180.                The above components can  be transformed into Earth-fixed-
  2181.                geocentric coordinates with the aid of a rotation matrix R
  2182.                which  shall  be  defined,  subsequently,  when coordinate
  2183.                systems and transformations are treated.  At that time, an
  2184.                innovative algorithm  for computing altitude  will also be
  2185.                developed.  For now, it suffices to write
  2186.  
  2187.                    ┌─  ─┐      ┌─  ─┐
  2188.                    │ q1 │      │ qe │
  2189.                    │ q2 │ =  R │ qn │
  2190.                    │ q3 │      │ qu │
  2191.                    └─  ─┘      └─  ─┘
  2192.  
  2193.  
  2194.  
  2195.  
  2196.                Trajectory Maker                                Drag Model
  2197.                Algorithm Supplement                                  4-15
  2198.                ──────────────────────────────────────────────────────────
  2199.                
  2200.  
  2201.  
  2202.  
  2203.  
  2204.  
  2205.  
  2206.  
  2207.  
  2208.  
  2209.  
  2210.  
  2211.  
  2212.  
  2213.  
  2214.  
  2215.  
  2216.  
  2217.  
  2218.  
  2219.  
  2220.  
  2221.  
  2222.  
  2223.  
  2224.  
  2225.  
  2226.  
  2227.  
  2228.  
  2229.  
  2230.  
  2231.  
  2232.  
  2233.  
  2234.  
  2235.  
  2236.  
  2237.  
  2238.  
  2239.  
  2240.  
  2241.  
  2242.  
  2243.  
  2244.  
  2245.  
  2246.  
  2247.  
  2248.                Fig. 4-8.  East-North Wind Velocity Components, Example 1.
  2249.  
  2250.  
  2251.  
  2252.  
  2253.                Trajectory Maker                                Drag Model
  2254.                Algorithm Supplement                                  4-16
  2255.                ──────────────────────────────────────────────────────────
  2256.                
  2257.  
  2258.  
  2259.  
  2260.  
  2261.  
  2262.  
  2263.  
  2264.  
  2265.  
  2266.  
  2267.  
  2268.  
  2269.  
  2270.  
  2271.  
  2272.  
  2273.  
  2274.  
  2275.  
  2276.  
  2277.  
  2278.  
  2279.  
  2280.  
  2281.  
  2282.  
  2283.  
  2284.  
  2285.  
  2286.  
  2287.  
  2288.  
  2289.  
  2290.  
  2291.  
  2292.  
  2293.  
  2294.  
  2295.  
  2296.  
  2297.  
  2298.  
  2299.  
  2300.  
  2301.  
  2302.  
  2303.  
  2304.  
  2305.                Fig. 4-9.  East-North Wind Velocity Components, Example 2.
  2306.  
  2307.  
  2308.  
  2309.  
  2310.                Trajectory Maker                                Drag Model
  2311.                Algorithm Supplement                                  4-17
  2312.                ──────────────────────────────────────────────────────────
  2313.                
  2314.  
  2315.  
  2316.                                        References
  2317.  
  2318.                1.  Menzner, et. al., Defining Constants, Equations,
  2319.                    and Abbreviated Tables of the 1975 U. S. Standard
  2320.                    Atmosphere, National Aeronautics and Space
  2321.                    Administration, Technical Report NASA TR R-459,
  2322.                    Washington D. C., May 1976.
  2323.  
  2324.                2.  Reynolds, H. B., Atmospheric Drag Effects on a Space
  2325.                    Based Radar Antenna, TRW Technical Memorandum,
  2326.                    January, 1991.
  2327.  
  2328.                3.  Thompson, R. F., Spline Interpolation on  a Digital
  2329.                    Computer, Goddard Space Flight Center, Technical
  2330.                    Report X-692-70-261, Greenbelt, Md., July 1970.
  2331.  
  2332.                4.  Waison, G. R., Prediction of Orbital Decay Rates Due
  2333.                    to Atmospheric Drag, TRW Program Technical Report,
  2334.                    3150-6020-RO-000, January 7, 1966.
  2335.  
  2336.                5.  Wertz, J. R., ed., Spacecraft Attitude Determination
  2337.                    and Control, D. Reidel Publishing Company, Boston,
  2338.                    1985.
  2339.  
  2340.  
  2341.  
  2342.  
  2343.  
  2344.                Trajectory Maker                        Integration Scheme
  2345.                Algorithm Supplement                                  5-1
  2346.                ──────────────────────────────────────────────────────────
  2347.                
  2348.  
  2349.  
  2350.                5.  NUMERICAL INTEGRATION SCHEME
  2351.  
  2352.                Substituting  the   atmospheric  drag   and  gravitational
  2353.                perturbation acceleration  components, into  the component
  2354.                form of the  first order vector  differential equations of
  2355.                motion, results  in the  following explicit  system of six
  2356.                first-order ordinary differential equations:
  2357.  
  2358.                    dx   .
  2359.                    ── = x
  2360.                    dt
  2361.  
  2362.                    dy   .
  2363.                    ── = y
  2364.                    dt
  2365.  
  2366.                    dz   .
  2367.                    ── = z
  2368.                    dt
  2369.                     .
  2370.                    dx     µx
  2371.                    ── = - ── + δGx + Dx
  2372.                    dt     r3
  2373.                     .
  2374.                    dy     µy
  2375.                    ── = - ── + δGy + Dy
  2376.                    dt     r3
  2377.                     .
  2378.                    dz     µz
  2379.                    ── = - ── + δGz + Dz
  2380.                    dt     r3
  2381.  
  2382.                If  the  perturbation  components  are  set  to  zero, the
  2383.                resulting equations are clearly recognizable as describing
  2384.                classical two-body motion.
  2385.  
  2386.                By defining the system state as the vector
  2387.                                . . .
  2388.                    x = [ x y z x y z ]T
  2389.  
  2390.                where the superscript T  denotes matrix transposition, the
  2391.                differential  equations   may  now   be  brought   into  a
  2392.                computational convenient vector form
  2393.  
  2394.                    dx
  2395.                    ── = f(x)
  2396.                    dt
  2397.  
  2398.  
  2399.  
  2400.  
  2401.                Trajectory Maker                        Integration Scheme
  2402.                Algorithm Supplement                                  5-2
  2403.                ──────────────────────────────────────────────────────────
  2404.                
  2405.  
  2406.  
  2407.                where f represents the vector-valued function of the state
  2408.                vector defined by the  right-hand-side of the differential
  2409.                equations.
  2410.  
  2411.                The differential  equations of motion  comprise an initial
  2412.                value system,  which may be  written in  more generally in
  2413.                the vector form including time,
  2414.  
  2415.                    dx
  2416.                    ── = f(t, x),   x(t0) = x0
  2417.                    dt
  2418.                                                               .  .  .
  2419.                where x = (x1, x2, x3, x4, x5, x6) ≡ (x, y, z, x, y, z).
  2420.  
  2421.                A numerical  scheme for integrating  initial value systems
  2422.                that  is  of  very  high-precision,  generally  stable and
  2423.                relatively  easy  to implement  is  Shanks'  explicit 8-12
  2424.                formulation  for  solutions to  differential  equations by
  2425.                evaluations of functions [1].
  2426.  
  2427.                This method is based on the classical Runge-Kutta formulas
  2428.                found in most textbooks on numerical analysis.
  2429.  
  2430.                The 8-12 means that the order of the integration scheme is
  2431.                eight (that  is, the truncation  error is  order nine) and
  2432.                twelve  evaluations (that  is,  stages) of  the right-hand
  2433.                side of the differential equations are performed.
  2434.  
  2435.                This scheme is a fixed  step integrator which is preferred
  2436.                for  trajectory  animation   graphics.   For  example,  an
  2437.                orbiting object may move farther  in its present time step
  2438.                than its preceding  time step. To  an observer, the object
  2439.                appears to move faster.
  2440.  
  2441.                An  advantage of  employing  a high-order  scheme  is that
  2442.                accuracy  may be  maintained  in considerable  fewer steps
  2443.                than the classical formulations.
  2444.  
  2445.                Given the state x = x(t) of the system at time t, then the
  2446.                state of the system at t + δt is given by
  2447.  
  2448.                    x(t + δt) = x(t) + δx
  2449.  
  2450.                where
  2451.  
  2452.                    δx ≈ ( 41 k1 + 216 k6 + 272 k7 + 27 k8 +  27 k9
  2453.  
  2454.                         + 36 k10 + 180 k11 + 41 k12 ) / 840 + O(δt9)
  2455.  
  2456.  
  2457.  
  2458.  
  2459.                Trajectory Maker                        Integration Scheme
  2460.                Algorithm Supplement                                  5-3
  2461.                ──────────────────────────────────────────────────────────
  2462.                
  2463.  
  2464.  
  2465.                Using the Einstein  summation convention with  σ = 1, ...,
  2466.                12 and µ ≤ σ, the coefficients may be defined by
  2467.  
  2468.                    kσ = f(t + tσ δt, x + sσµ kµ ) δt
  2469.  
  2470.                The  non-zero  constants  tσ  and  sσµ  appearing  in this
  2471.                expression are defined in Table 5-1 below.
  2472.  
  2473.  
  2474.                 Table 5-1. Shanks' 8-12 Numerical Integration Constants.
  2475.  
  2476.                ╒════════════════════════════════════════════════════════╕
  2477.                │ t2 = 1/9     t3 = 1/6   t4 = 1/4    t5 = 1/10          │
  2478.                │ t6 = 1/6     t7 = 1/2   t8 = 2/3    t9 = 1/3           │
  2479.                │ t10= 5/6     t11= 5/6   t12= 1                         │
  2480.                ├────────────────────────────────────────────────────────┤
  2481.                │ s11 = 1/9                                              │
  2482.                │ s21 = 1/24         s22 = 1/8                           │
  2483.                │ s31 = 1/16         s33 = 3/16                          │
  2484.                │ s41 = 29/500       s43 = 33/500      s44 = -3/125      │
  2485.                │ s51 = 11/324       s54 = 1/243       s55 = 125/972     │
  2486.                ├────────────────────────────────────────────────────────┤
  2487.                │ s61 = -7/12        s64 = 19/9        s65 = 125/36      │
  2488.                │ s66 = -9/2                                             │
  2489.                ├────────────────────────────────────────────────────────┤
  2490.                │ s71 = -10/81       s74 = -32/243     s75 = 125/243     │
  2491.                │ s77 = 11/27                                            │
  2492.                ├────────────────────────────────────────────────────────┤
  2493.                │ s81 = 1175/324     s84 = -32/3        s85 = -3125/162  │
  2494.                │ s86 = 26           s87 = 121/162      s88 = -1/12      │
  2495.                ├────────────────────────────────────────────────────────┤
  2496.                │ s91 = 293/324      s94 = -71/27       s95 = -1375/324  │
  2497.                │ s96 = 51/9         s97 = -59/162      s98 = 1/2        │
  2498.                │ s99 = 1                                                │
  2499.                ├────────────────────────────────────────────────────────┤
  2500.                │ s101 = 1303/1620   s104 = -71/27      s105 = -1375/324 │
  2501.                │ s106 = 37/6        s107 = 103/162     s1010= 1/10      │
  2502.                ├────────────────────────────────────────────────────────┤
  2503.                │ s111 = -955/492    s114 = 2560/369    s115 = 8125/738  │
  2504.                │ s116 = -612/41     s117 = 7/82        s118 = -27/164   │
  2505.                │ s119 = -18/41      s1110= -12/41      s1111=  30/41    │
  2506.                └────────────────────────────────────────────────────────┘
  2507.  
  2508.  
  2509.                Then, given the initial state x0 of the system at time t0,
  2510.                the time increment δt and  the final time tf, the solution
  2511.                may  be  generated  (or  propagated)  with  the  following
  2512.                recursion scheme:
  2513.  
  2514.  
  2515.  
  2516.  
  2517.                Trajectory Maker                        Integration Scheme
  2518.                Algorithm Supplement                                  5-4
  2519.                ──────────────────────────────────────────────────────────
  2520.                
  2521.  
  2522.  
  2523.                    x(tm) = x(t0 + m δt);   m = 0, 1, ..., (tf - t0)/δt
  2524.  
  2525.                A  numerical  scheme   that  integrates  the  differential
  2526.                equations of motion  directly, such as  this one, is known
  2527.                as a Cowell integration scheme.  The scheme doesn't really
  2528.                care  what kind  of  orbit it  is  generating, and  it can
  2529.                handle   deviations   from   classical   two-body  motion.
  2530.                Moreover,  the user  has  complete control  concerning its
  2531.                accuracy with his choice of propagation step-size.
  2532.  
  2533.                Figure  5-1  illustrates  the   common  logarithm  of  the
  2534.                relative  error  in  Shanks'  8-12  numerical  integration
  2535.                scheme.
  2536.  
  2537.                The error is based on a 200 nm unperturbed orbit.  In this
  2538.                instance, the  Earth is  modeled as  a sphere  so that, at
  2539.                least theoretically,  the altitude  is a  constant 200 nm.
  2540.                The  computed altitude  will, of  course, differ  from the
  2541.                true altitude  depending on the  integration step-size and
  2542.                the number of integration steps.
  2543.  
  2544.                The  figure illustrates  the error  as  a function  of the
  2545.                final time  in days of  86400 seconds.  (Equivalent to the
  2546.                duration of the orbit since the initial time was zero.)
  2547.  
  2548.                The relative error is shown because of its connection with
  2549.                the number  of significant  figures and  independence from
  2550.                any particular units  employed.  From elementary numerical
  2551.                analysis, it  is known that  if the relative  error in any
  2552.                number  is not  greater  than 0.5  X  10-N, the  number is
  2553.                certainly correct to N significant figures.
  2554.  
  2555.                It is  seen, for example,  that if an  orbit is propagated
  2556.                with a step size of 300 seconds (5 minutes) for a duration
  2557.                of  7 days  (2016 steps),  the  logarithm of  the relative
  2558.                error is about -5.5.  Since
  2559.  
  2560.                    log(0.5 X 10-N) ≈ -(0.301 + N)
  2561.  
  2562.                and  -5.5  <  -5.301, the  computation  is  accurate  to 5
  2563.                significant figures, which is  equivalent to about 0.01 nm
  2564.                for a 200.00 nm orbit.
  2565.  
  2566.                It is emphasized that the results obtained from the figure
  2567.                are  based on  an unperturbed  orbit.  However,  they will
  2568.                generally be valid for gravity perturbed orbits as well.
  2569.  
  2570.  
  2571.  
  2572.  
  2573.                Trajectory Maker                        Integration Scheme
  2574.                Algorithm Supplement                                  5-5
  2575.                ──────────────────────────────────────────────────────────
  2576.                
  2577.  
  2578.  
  2579.                Things  are  quite  different  for  objects  entering  the
  2580.                Earth's atmosphere.  Such an object  experiences a violent
  2581.                transient in  its acceleration  profile that  can tax most
  2582.                any integrator's numerical integrity.
  2583.  
  2584.                Several factors are  involved. The speed  of the object is
  2585.                squared in the drag force  equation, the density varies by
  2586.                orders of  magnitude becoming large  in the  100 km region
  2587.                and the drag force is sensitive to the ballistic coeffi-
  2588.                cient. (Small ballistic coefficients result in larger drag
  2589.                forces than large ballistic coefficients.)
  2590.  
  2591.  
  2592.  
  2593.  
  2594.                Trajectory Maker                        Integration Scheme
  2595.                Algorithm Supplement                                  5-6
  2596.                ──────────────────────────────────────────────────────────
  2597.                
  2598.  
  2599.  
  2600.  
  2601.  
  2602.  
  2603.  
  2604.  
  2605.  
  2606.  
  2607.  
  2608.  
  2609.  
  2610.  
  2611.  
  2612.  
  2613.  
  2614.  
  2615.  
  2616.  
  2617.  
  2618.  
  2619.  
  2620.  
  2621.  
  2622.  
  2623.  
  2624.  
  2625.  
  2626.  
  2627.  
  2628.  
  2629.  
  2630.  
  2631.  
  2632.  
  2633.  
  2634.  
  2635.  
  2636.  
  2637.  
  2638.  
  2639.  
  2640.  
  2641.  
  2642.  
  2643.  
  2644.  
  2645.  
  2646.                Figure 5-1.  Numerical Integration Error
  2647.  
  2648.  
  2649.  
  2650.  
  2651.                Trajectory Maker                        Integration Scheme
  2652.                Algorithm Supplement                                  5-7
  2653.                ──────────────────────────────────────────────────────────
  2654.                
  2655.  
  2656.  
  2657.                                        References
  2658.  
  2659.                1.  Shanks, B. F., Solutions of Differential Equations by
  2660.                    Evaluations of Functions, Math. Comp. 20 (1966) pp.
  2661.                    21-38.
  2662.  
  2663.  
  2664.  
  2665.  
  2666.                Trajectory Maker                        Coordinate Systems
  2667.                Algorithm Supplement                                  6-1
  2668.                ──────────────────────────────────────────────────────────
  2669.                
  2670.  
  2671.  
  2672.                6.  COORDINATE SYSTEMS AND TRANSFORMATIONS
  2673.  
  2674.                We have seen  that in order to  propagate the equations of
  2675.                motion  with  a  Cowell integration  scheme,  all  that is
  2676.                required  is  an  initial state  vector  given  in  Earth-
  2677.                centered-inertial coordinates,
  2678.                             .  .  .
  2679.                    [x, y, z, x, y, z]T
  2680.  
  2681.                at some epoch  t0.  Rectangular coordinates  such as these
  2682.                are  difficult  to visualize.   Nevertheless,  the program
  2683.                should  include the  option for   specifying input  in the
  2684.                form of ECI coordinates, because these coordinates will be
  2685.                known for  many applications; in  particular for ballistic
  2686.                missile trajectories.
  2687.  
  2688.                However, there  exist other coordinate  systems which will
  2689.                simplify  our  visualization  process  when  it  comes  to
  2690.                describing  an  initial trajectory  with  input  data, and
  2691.                interpreting various characteristics with output data.
  2692.  
  2693.  
  2694.                6.1  Time Considerations
  2695.  
  2696.                A  fundamental  parameter  that  is  needed  to  define an
  2697.                inertial reference system  is the Greenwich  hour angle of
  2698.                the  vernal equinox  at  epoch or  orbit  insertion.  This
  2699.                angle is  equivalent to  Greenwich sidereal  time at epoch
  2700.                and serves to fix a satellite ground trace relative to the
  2701.                rotating earth.
  2702.  
  2703.                We suppose that  the initial time,  or epoch, is specified
  2704.                by the year,  I, month, J,  day, K, and  universal time in
  2705.                hours, minutes, and seconds  (past midnight).  These units
  2706.                may  be converted  to the  Julian  date JD,  measured from
  2707.                January 0.0 by the clever expression in Wertz [4],
  2708.  
  2709.                 JD = K - 32075 + 1461*( I + 4800 + (J - 14)/12)/4
  2710.                                 + 367*( J -    2 - (J - 14)/12*12)/12
  2711.                                   - 3*((I + 4900 + (J - 14)/12)/100)/4
  2712.                                   - 0.50
  2713.  
  2714.                where  I,  J,  K  are  defined  to  be  integers  and  the
  2715.                arithmetic  is  performed  in   integer  mode.   That  is,
  2716.                floating point  numbers are truncated  at each  stage of a
  2717.                computation, such as done by a computer language such as C
  2718.                or FORTRAN.
  2719.  
  2720.  
  2721.  
  2722.  
  2723.                Trajectory Maker                        Coordinate Systems
  2724.                Algorithm Supplement                                  6-2
  2725.                ──────────────────────────────────────────────────────────
  2726.                
  2727.  
  2728.  
  2729.                The right ascension of the Sun is obtained using Newcomb's
  2730.                equation [1]
  2731.  
  2732.                    Ru = 279.6910 + 36000.7689 T + 0.0004 T2   (degrees)
  2733.  
  2734.                where T  is the number  of Julian centuries  of 36525 days
  2735.                elapsed since noon Greenwich mean time on January 0, 1900.
  2736.                This value is computed from
  2737.  
  2738.                    T = (JD - JD0)/36525
  2739.  
  2740.                where JD0  is the Julian  day from Greenwich  mean noon on
  2741.                January 0.5, 1900 which is 2,415,020.0.
  2742.  
  2743.                Universal time in degrees is
  2744.  
  2745.                    UT = 15*(hours + minutes/60 + seconds/3600)
  2746.  
  2747.                So that the Greenwich sidereal time is
  2748.  
  2749.                    GST = Ru + UT - 180   (degrees)
  2750.  
  2751.  
  2752.                6.2  Orbital Elements
  2753.  
  2754.                A  convenient method  for defining  an  orbit is  with its
  2755.                orbital elements.  These elements describe its size, shape
  2756.                and orientation.   They are:  semimajor  axis a, eccentri-
  2757.                city e,  inclination i, longitude of the ascending node Ω,
  2758.                argument of  perigee w,  and some  time when  the orbiting
  2759.                object passes perigee τ.
  2760.  
  2761.                This latter element is also  equivalent to mean anomaly M,
  2762.                but  when  specifying non-elliptical  orbits,  it  is more
  2763.                convenient to specify  time past perigee  and compute mean
  2764.                anomaly from the formula
  2765.  
  2766.                    M = nτ
  2767.  
  2768.                where n is the mean motion given by
  2769.  
  2770.                    n = (µ a3)1/2
  2771.  
  2772.                This parameter  is the  same as  2π/P for  elliptic orbits
  2773.                having period P.
  2774.  
  2775.                It is a  known from two-body  theory that orbital elements
  2776.                at an  instant in  time, in  particular at  the epoch, are
  2777.  
  2778.  
  2779.  
  2780.  
  2781.                Trajectory Maker                        Coordinate Systems
  2782.                Algorithm Supplement                                  6-3
  2783.                ──────────────────────────────────────────────────────────
  2784.                
  2785.  
  2786.  
  2787.                integration  constants  of the  differential  equations of
  2788.                motion,  and   these  constants  are   equivalent  to  the
  2789.                instantaneous   rectangular   coordinates   by   a   given
  2790.                coordinate  transformation.  This  transformation  is from
  2791.                Roy [3].
  2792.  
  2793.                A geocentric, inertial coordinate system is defined in the
  2794.                orbit plane  with principle axis  x* along  the major axis
  2795.                positive  in the  direction  of perigee,  another  axis z*
  2796.                normal to  the orbit  plane in  the direction  of positive
  2797.                angular momentum and  a third axis y*  defined so that the
  2798.                system is right handed.
  2799.  
  2800.                From the orientation angles Ω, w, i, the direction cosines
  2801.                of the  x*, y*,  z* two  planar axes  relative to  the ECI
  2802.                system may be obtained.  They are
  2803.  
  2804.                    l1 = cos Ω cos w - sin Ω sin w cos i
  2805.                    m1 = sin Ω cos w + cos Ω sin w cos i
  2806.                    n1 = sin w sin i
  2807.  
  2808.                    l2 = -cos Ω cos w - sin Ω sin w cos i
  2809.                    m2 = -sin Ω cos w + cos Ω sin w cos i
  2810.                    n2 = cos w sin i
  2811.  
  2812.                These  cosines are  independent  of the  specific  type of
  2813.                orbit be it elliptic, hyperbolic or parabolic.
  2814.  
  2815.  
  2816.                Elliptical Case
  2817.                ───────────────
  2818.                
  2819.                The ECI state vector for the elliptical case is given by
  2820.  
  2821.                    x = al1 cos E + bl2 sin E - ael1
  2822.  
  2823.                    y = am1 cos E + bm2 sin E - aem1
  2824.  
  2825.                    z = an1 cos E + bn2 sin E - aen1
  2826.                    .
  2827.                    x = (bl2 cos E - al1 sin E)(na/r)
  2828.                    .
  2829.                    y = (bm2 cos E - am1 sin E)(na/r)
  2830.                    .
  2831.                    z = (bn2 cos E - an1 sin E)(na/r)
  2832.  
  2833.                where b is the orbit's semiminor axis given by
  2834.  
  2835.  
  2836.  
  2837.  
  2838.                Trajectory Maker                        Coordinate Systems
  2839.                Algorithm Supplement                                  6-4
  2840.                ──────────────────────────────────────────────────────────
  2841.                
  2842.  
  2843.  
  2844.                    b = a(1 - e2)1/2
  2845.  
  2846.                and r is objects' radius vector magnitude given by
  2847.  
  2848.                    r = a(1 - e cos E)
  2849.  
  2850.                The  parameter   E  is  the   object's  eccentric  anomaly
  2851.                obtainable from Kepler equation
  2852.  
  2853.                    E - e sin E = M
  2854.  
  2855.                by employing Newton's method for finding roots.
  2856.  
  2857.  
  2858.                Hyperbolic Case
  2859.                ───────────────
  2860.                
  2861.                The ECI state vector for the hyperbolic case is given by
  2862.  
  2863.                    x = -al1 cosh F + bl2 sinh F + ael1
  2864.  
  2865.                    y = -am1 cosh F + bm2 sinh F + aem1
  2866.  
  2867.                    z = -an1 cosh F + bn2 sinh F + aen1
  2868.                    .
  2869.                    x = (bl2 cosh F - al1 sinh F)(na/r)
  2870.                    .
  2871.                    y = (bm2 cosh F - am1 sinh F)(na/r)
  2872.                    .
  2873.                    z = (bn2 cosh F - an1 sinh F)(na/r)
  2874.  
  2875.  
  2876.                where b is the orbit's semiminor axis given by
  2877.  
  2878.                    b = a(e2 - 1)1/2
  2879.  
  2880.                and r is objects' radius vector magnitude given by
  2881.  
  2882.                    r = a(F cosh F - 1)
  2883.  
  2884.                The parameter F is the objects' hyperbolic areal parameter
  2885.                obtainable from the hyperbolic analog of Kepler's equation
  2886.  
  2887.                    F sinh F - F = M
  2888.  
  2889.                also, by employing Newton's method for finding roots.
  2890.  
  2891.  
  2892.  
  2893.  
  2894.                Trajectory Maker                        Coordinate Systems
  2895.                Algorithm Supplement                                  6-5
  2896.                ──────────────────────────────────────────────────────────
  2897.                
  2898.  
  2899.  
  2900.                Parabolic Case
  2901.                ───────────────
  2902.  
  2903.                The ECI state vector for the parabolic case is given by
  2904.  
  2905.                    x = l1r x* + l2r y*
  2906.  
  2907.                    y = m1r x* + m2r y*
  2908.  
  2909.                    z = n1r x* + n2r y*
  2910.                    .
  2911.                    x = l1 x* + l2 y*
  2912.                    .
  2913.                    y = m1 x* + m2 y*
  2914.                    .
  2915.                    z = n1 x* + n2 y*
  2916.  
  2917.                where
  2918.  
  2919.                    x* = √(µp) ( sin f* cos f* /p - sin f* /r)
  2920.  
  2921.                    y* = √(µp) ( sin f* sin f* /p - cos f* /r)
  2922.  
  2923.                For the parabolic  case, the semimajor  axis a is defined,
  2924.                herein,  as the  object's perigee  distance and  the conic
  2925.                parameter is computed from
  2926.  
  2927.                    p = 2a
  2928.  
  2929.                The parameter f*  is the object's true  anomaly and may be
  2930.                found by solving Barker's cubic equation
  2931.                                  _
  2932.                    D + 1/3 D3 = 2n(t - τ)
  2933.  
  2934.                where
  2935.  
  2936.                    D = tan (f*/2)
  2937.  
  2938.                and
  2939.                    _
  2940.                    n2 p3 = µ
  2941.  
  2942.                Details of the solution are given in Roy [3].
  2943.  
  2944.  
  2945.  
  2946.  
  2947.                Trajectory Maker                        Coordinate Systems
  2948.                Algorithm Supplement                                  6-6
  2949.                ──────────────────────────────────────────────────────────
  2950.                
  2951.  
  2952.  
  2953.                6.3  Geospherical Inertial Coordinates
  2954.  
  2955.                Geospherical    inertial    coordinates    are   spherical
  2956.                coordinates  having their  origin  at the  Earth's center.
  2957.                These  coordinates consist  of an  objects range  r, right
  2958.                ascension α and declination  δ.  Declination is equivalent
  2959.                to geocentric latitude, until now, denoted by φ'.
  2960.  
  2961.                These coordinates are defined by
  2962.  
  2963.                    r = √(x2 + y2 + z2)
  2964.  
  2965.                    δ = arcsin(z/r)
  2966.  
  2967.                    α = arctan(y/x)
  2968.  
  2969.                In addition, the following flight parameters are computed:
  2970.  
  2971.                Inertial velocity
  2972.  
  2973.                          .    .    .
  2974.                    v = √(x2 + y2 + z2)
  2975.  
  2976.                Flight path angle
  2977.  
  2978.                    Γ = arcsin(x'/v)
  2979.  
  2980.                Velocity azimuth angle
  2981.  
  2982.                    σ = arctan(y'/z')
  2983.  
  2984.                where
  2985.  
  2986.                 ┌─  ─┐   ┌─                                   ─┐ ┌─ ─┐
  2987.                 │ x' │   │  cos δ cos α    cos δ sin α   sin δ │ │ x │
  2988.                 │    │   │                                     │ │   │
  2989.                 │ y' │ = │    -sin α          cos α        0   │ │ y │
  2990.                 │    │   │                                     │ │   │
  2991.                 │ z' │   │ -sin δ cos α   -sin δ sin α   cos δ │ │ z │
  2992.                 └─  ─┘   └─                                   ─┘ └─ ─┘
  2993.  
  2994.                All the remaining coordinates used for the software system
  2995.                are relative to the rotating Earth.
  2996.  
  2997.  
  2998.                6.4  Earth-Fixed-Geocentric Coordinates
  2999.  
  3000.                Earth-fixed-geocentric  (EFG) coordinates  have  origin at
  3001.  
  3002.  
  3003.  
  3004.  
  3005.                Trajectory Maker                        Coordinate Systems
  3006.                Algorithm Supplement                                  6-7
  3007.                ──────────────────────────────────────────────────────────
  3008.                
  3009.  
  3010.  
  3011.                the Earth's  geocenter, fundamental  plane coincident with
  3012.                the equatorial plane  and principle axis  in the direction
  3013.                of the Greenwich meridian.   This coordinate system, being
  3014.                fixed  relative  to  the  Earth  is  actually  a  rotating
  3015.                coordinate system.
  3016.  
  3017.                Letting e, f, g denote  these coordinates, the e-axis lies
  3018.                in the equatorial  plane and is  positive in the direction
  3019.                of Greenwich;  the g-axis  is coincident  with the Earth's
  3020.                polar axis, i.e.,  its spin axis; and  the f-axis lying in
  3021.                the equatorial plane such that the system is right-handed.
  3022.  
  3023.                This coordinate system is related to the inertial frame by
  3024.                the transformation
  3025.  
  3026.                    e =  x cos(gst) + y sin(gst)
  3027.  
  3028.                    f = -x sin(gst) + y cos(gst)
  3029.  
  3030.                    g =  z
  3031.  
  3032.                where gst is the current  Greenwich sidereal time.  If gha
  3033.                is the Greenwich sidereal time at epoch, then
  3034.  
  3035.                    gst= gha + Ωt  (mod 360)
  3036.  
  3037.                where t is the time past epoch.
  3038.  
  3039.  
  3040.                6.5  Local Topocentric Coordinates
  3041.  
  3042.                Local topocentric  coordinates (UVW)  have origin  at some
  3043.                local site  on the surface  of the  Earth, and fundamental
  3044.                plane tangent to the surface  at the site.  This reference
  3045.                system is also a rotating coordinate system.
  3046.  
  3047.                Letting u, v, w denote these coordinates, where the w-axis
  3048.                is  normal  to the  tangent  plane at  the  site, positive
  3049.                upwards;  the  principle  u-axis is  lies  in  the tangent
  3050.                plane, positive at  in an azimuthal  direction   ; and the
  3051.                v-axis lies in  the tangent plane such  that the system is
  3052.                right-handed.
  3053.  
  3054.                It is shown in O'Connor  [2], that the transformation from
  3055.                earth-fixed-geocentric  coordinates  to  local topocentric
  3056.                coordinates is given by the transformation
  3057.  
  3058.  
  3059.  
  3060.  
  3061.                Trajectory Maker                        Coordinate Systems
  3062.                Algorithm Supplement                                  6-8
  3063.                ──────────────────────────────────────────────────────────
  3064.                
  3065.  
  3066.  
  3067.                    ┌─ ─┐       ┌─      ─┐
  3068.                    │ u │       │ e - e0 │
  3069.                    │ v │ =  ABC│ f - f0 │
  3070.                    │ w │       │ g - g0 │
  3071.                    └─ ─┘       └─      ─┘
  3072.  
  3073.                where the rotation matrices are
  3074.  
  3075.                          ┌─                    ─┐
  3076.                          │  sin      cos      0 │
  3077.                    A =   │ -cos      sin      0 │
  3078.                          │    0        0      1 │
  3079.                          └─                    ─┘
  3080.  
  3081.                          ┌─                    ─┐
  3082.                          │ 1      0        0    │
  3083.                    B =   │ 0    sin φ    cos φ  │
  3084.                          │ 0   -cos φ    sin φ  │
  3085.                          └─                    ─┘
  3086.  
  3087.                          ┌─                    ─┐
  3088.                          │ -sin      cos      0 │
  3089.                    C =   │ -cos     -sin      0 │
  3090.                          │    0        0      1 │
  3091.                          └─                    ─┘
  3092.  
  3093.                where φ,   denote the latitude  and longitude of the site,
  3094.                respectively, and e0, f0, g0, are its EFG coordinates.
  3095.  
  3096.                These latter coordinates may be computed directly from
  3097.  
  3098.                         ┌─                     ─┐
  3099.                         │        a              │
  3100.                    e0 = │ ────────────────  + h │ cos φ cos *
  3101.                         │ √(1 - e2 sin2 φ)      │
  3102.                         └─                     ─┘
  3103.  
  3104.                         ┌─                     ─┐
  3105.                         │        a              │
  3106.                    f0 = │ ────────────────  + h │ cos φ sin *
  3107.                         │ √(1 - e2 sin2 φ)      │
  3108.                         └─                     ─┘
  3109.  
  3110.                         ┌─                     ─┐
  3111.                         │    a(1 - e2)          │
  3112.                    g0 = │ ────────────────  + h │ sin φ
  3113.                         │ √(1 - e2 sin2 φ)      │
  3114.                         └─                     ─┘
  3115.  
  3116.  
  3117.  
  3118.  
  3119.                Trajectory Maker                        Coordinate Systems
  3120.                Algorithm Supplement                                  6-9
  3121.                ──────────────────────────────────────────────────────────
  3122.                
  3123.  
  3124.  
  3125.                If the azimuth angle is 90°, this system is referred to as
  3126.                the  East-North-Up coordinate  system.   In this  case the
  3127.                azimuthal matrix is the unit matrix
  3128.  
  3129.                         ┌─           ─┐
  3130.                         │ 1    0    0 │
  3131.                    A =  │ 0    1    0 │
  3132.                         │ 0    0    1 │
  3133.                         └─           ─┘
  3134.  
  3135.                and the transformation becomes
  3136.  
  3137.                    ┌─  ─┐      ┌─      ─┐
  3138.                    │ ue │      │ e - e0 │
  3139.                    │ vn │ =  BC│ f - f0 │
  3140.                    │ wu │      │ g - g0 │
  3141.                    └─  ─┘      └─      ─┘
  3142.  
  3143.                where
  3144.  
  3145.                          ┌─                                      ─┐
  3146.                          │        -sin            cos        0    │
  3147.                    BC=   │  -sin φ cos     -sin φ sin        0    │
  3148.                          │   cos φ cos      cos φ sin      sin φ  │
  3149.                          └─                                      ─┘
  3150.  
  3151.                By differentiating  and solving for  the e,  f, g velocity
  3152.                components, we find that
  3153.  
  3154.                     .  .  .                .  .  .
  3155.                   [ e  f  g ]T =   (BC)T [ u  v  w ]T
  3156.  
  3157.                These  components   may  be   identified  with   the  wind
  3158.                components in Section 4.  That is  R= (BC)T so that
  3159.  
  3160.                   ┌   ─┐      ┌─  ─┐
  3161.                   │ q1 │      │ qe │
  3162.                   │ q2 │ =  R │ qn │
  3163.                   │ q3 │      │ qu │
  3164.                   └   ─┘      └─  ─┘
  3165.  
  3166.                Carrying out the matrix operations, we have
  3167.  
  3168.                    q1 = -[qe sin   + qn cos * sin φ]
  3169.  
  3170.                    q2 =  [qe cos   - qn sin * sin φ]
  3171.  
  3172.                    q3 =  qn cos φ
  3173.  
  3174.  
  3175.  
  3176.  
  3177.                Trajectory Maker                        Coordinate Systems
  3178.                Algorithm Supplement                                  6-10
  3179.                ──────────────────────────────────────────────────────────
  3180.                
  3181.  
  3182.  
  3183.                which are  the inertial wind  velocity components required
  3184.                for the equations of motion.
  3185.  
  3186.  
  3187.                6.6  Common Radar Coordinates
  3188.  
  3189.                Radar coordinates  (RAE) are  used for  measurements of at
  3190.                some particular geographic location or site.  They consist
  3191.                of an object's  range, azimuth and  elevation.  Range R is
  3192.                the distance from the site to the object; azimuth A is the
  3193.                angle between  the object  and the  radar site's meridian,
  3194.                measured  clockwise  from   the  northern  direction;  and
  3195.                elevation  E  is the  angle  above or  below  the horizon.
  3196.                These coordinates are given by the following expressions
  3197.  
  3198.                    R = √(ue2 +vn2 + wu2)
  3199.  
  3200.                    A = arctan(ue/vn)
  3201.  
  3202.                    E = arcsin(wn/R)
  3203.  
  3204.                where the subscripts e, n, u denote the east, north and up
  3205.                components, respectively.
  3206.  
  3207.                In  addition  range rate  is  obtained  by differentiating
  3208.                range.  We have,
  3209.  
  3210.                    .      .      .       .
  3211.                    R = (ueue + vnvn +  wuwu)/R
  3212.  
  3213.  
  3214.                6.7  Geodetic Coordinates and Computation of Altitude
  3215.  
  3216.                Altitude of  an orbiting  object is  its height  above the
  3217.                Earth's surface.  If the Earth  were a sphere of radius of
  3218.                radius a, then the  altitude of an object  at a distance r
  3219.                from the geocenter would be
  3220.  
  3221.                    h = r - a
  3222.  
  3223.                But the equatorial and polar radii differ by nearly 12 nm;
  3224.                The  polar  radius  being  the  smaller  of  the  two.   A
  3225.                parameter  commonly used  for measuring  the shape  of the
  3226.                Earth is its flattening defined by
  3227.  
  3228.                    f = (a - b)/a
  3229.  
  3230.                where,  in  this  section, a   and  b  denote  the Earth's
  3231.  
  3232.  
  3233.  
  3234.  
  3235.                Trajectory Maker                        Coordinate Systems
  3236.                Algorithm Supplement                                  6-11
  3237.                ──────────────────────────────────────────────────────────
  3238.                
  3239.  
  3240.  
  3241.                equatorial and  polar radii,  respectively.  Flattening is
  3242.                usually provided  as a  reciprocal, for  example the value
  3243.                used herein, is the WGS-72 value of 1/298.26.
  3244.  
  3245.                An improved altitude model is
  3246.  
  3247.                    h = r - a(1 - f sin φ')
  3248.  
  3249.                where,  as  before  φ',  denotes  the  declination  of the
  3250.                object.  The  declination   is  equivalent  to  geocentric
  3251.                latitude.
  3252.  
  3253.                Both  of  these  expressions are  used  by  the trajectory
  3254.                software, but  are principally  used for  estimates in the
  3255.                computer graphics.
  3256.  
  3257.                A great deal of effort was expended in order to accurately
  3258.                model the density profile, so we should expect an accurate
  3259.                altitude model from which density is obtained during orbit
  3260.                propagation.
  3261.  
  3262.                An altitude model that will  achieve any desired degree of
  3263.                precision that is also extremely fast; indeed the one used
  3264.                in this paper, is an unpublished algorithm conceived by L.
  3265.                H. Ivy.   A by-product of  his algorithm  are the trigono-
  3266.                metric functions  required to construct  the geocentric to
  3267.                topocentric rotation matrix, R sited in previous sections.
  3268.  
  3269.                Consider  an  instantaneous  cross-section  of  the  Earth
  3270.                defined by  the plane  containing its  polar axis  and the
  3271.                mass m.   Define a two-dimensional  coordinate system with
  3272.                u-axis formed by  the intersection of  this plane with the
  3273.                equatorial  plane, and  v-axis  coincident with  the polar
  3274.                axis.
  3275.  
  3276.                The cross-section is bounded  by an ellipse whose equation
  3277.                is
  3278.  
  3279.                    u2   v2
  3280.                    ── +  ── = 1
  3281.                    a2   b2
  3282.  
  3283.                Since the eccentricity of the ellipse is given by
  3284.  
  3285.                    b2 = a2(1 - e2)
  3286.  
  3287.                the its equation may be written
  3288.  
  3289.  
  3290.  
  3291.  
  3292.                Trajectory Maker                        Coordinate Systems
  3293.                Algorithm Supplement                                  6-12
  3294.                ──────────────────────────────────────────────────────────
  3295.                
  3296.  
  3297.  
  3298.                    u2(1 - e2) + v2 = a2(1 - e2)
  3299.  
  3300.                Now, the tangent of the geodetic  latitude φ of a point on
  3301.                the ellipse is the negative reciprocal of the slope of the
  3302.                line tangent to the ellipse  at that point.  It follows by
  3303.                differentiation that
  3304.  
  3305.                                1      v
  3306.                    tan φ =  ──────── ───
  3307.                             (1 - e2)  u
  3308.  
  3309.                Consequently,
  3310.  
  3311.                              a cos φ
  3312.                    u =  ────────────────
  3313.                         √(1 - e2 sin2 φ)
  3314.  
  3315.                          a(1 - e2) sin φ
  3316.                    v =  ────────────────
  3317.                         √(1 - e2 sin2 φ)
  3318.  
  3319.                These coordinates  are the coordinates  of a  point on the
  3320.                ellipse whose geodetic latitude is φ.
  3321.  
  3322.                On  the  other hand,  the  coordinates  of a  point  at an
  3323.                altitude h whose geodetic latitude is φ are
  3324.  
  3325.                    uh= u + h cos φ
  3326.  
  3327.                    vh= v + h sin φ
  3328.  
  3329.                The geocentric latitude φ' of this point is given by
  3330.  
  3331.                             ┌─                                ─┐
  3332.                             │                   e2             │
  3333.                    tan φ' = │ 1 -  ─────────────────────────── │ tan φ
  3334.                             │      1 + (h/a) √(1 - e2 sin2 φ ) │
  3335.                             └─                                ─┘
  3336.  
  3337.                One can also find the following expression for altitude
  3338.  
  3339.                    h= uh cos φ  + vh sin φ - a√(1 - e2 sin2 φ)
  3340.  
  3341.                These two equations form the basis of Ivy's algorithm.
  3342.  
  3343.                In essence,  the algorithm  is an  iterative process which
  3344.                uses these equations starting  with geocentric latitude as
  3345.                an estimate for geodetic latitude.
  3346.  
  3347.  
  3348.  
  3349.  
  3350.                Trajectory Maker                        Coordinate Systems
  3351.                Algorithm Supplement                                  6-13
  3352.                ──────────────────────────────────────────────────────────
  3353.                
  3354.  
  3355.  
  3356.                To complete the algorithm, Ivy  defines a miss distance as
  3357.                the  distance between  the  current estimate  for geodetic
  3358.                latitude φi and the line normal to the reference ellipsoid
  3359.                corresponding to the true geodetic coordinates.
  3360.  
  3361.                Let ui, vi be  a point on the  ellipsoid whose latitude is
  3362.                φi.  Then the equation of the line normal to the ellipsoid
  3363.                at this point, expressed in normal form, is
  3364.  
  3365.                   (v - vi)cos φi - (u - ui)sin φi = 0
  3366.  
  3367.                The distance from the line to the point u, v is
  3368.  
  3369.                    p = │ (u - ui)sin φi - (v - vi)cos φi │
  3370.  
  3371.                Using (3) and (4) with u = ui, v = vi, we have the
  3372.                convergence criterion,
  3373.  
  3374.  
  3375.                        │                          ae2sin φi cos φi │
  3376.                    p = │ uh sin φi - vh cos φi - ───────────────── │
  3377.                        │                          √(1 - e2 sin2 φ) │
  3378.  
  3379.  
  3380.                The algorithm can now be summarized.
  3381.  
  3382.  
  3383.                Algorithm Summary
  3384.                ─────────────────
  3385.                
  3386.                1. Input the coordinates of the object: x, y, z
  3387.  
  3388.                2. Compute:  √(x2 + y2)
  3389.  
  3390.                3. Compute:  tan φ' = z/√(x2 + y2)
  3391.  
  3392.                4. Replace:  tan φ = tan φ'
  3393.  
  3394.                5. Compute:  cos2 φ = 1/(1 + tan2 φ)
  3395.  
  3396.                6. Compute:  sin2 φ = 1 - cos2 φ
  3397.  
  3398.                7. Compute:  rdcl = √(1 - e2 sin2 φ)
  3399.  
  3400.                8. Compute:  cos φ = √cos2 φ
  3401.  
  3402.                9. Compute:  sin φ = cos φ tan φ
  3403.  
  3404.  
  3405.  
  3406.  
  3407.                Trajectory Maker                        Coordinate Systems
  3408.                Algorithm Supplement                                  6-14
  3409.                ──────────────────────────────────────────────────────────
  3410.                
  3411.  
  3412.  
  3413.                10. Compute: h = x cos φ + y sin φ - a rdcl
  3414.  
  3415.                11. Compute: p = │x sin φ - y cos φ - ae2sin φ cos φ/rdcl│
  3416.  
  3417.                12. Test: If p is less than some error criterion, end
  3418.                    algorithm.  Otherwise,
  3419.  
  3420.                13. Compute: D = 1 + (h/a) rdcl
  3421.  
  3422.                14. Compute: tan φ = tan φ '/(1 - e2/D)
  3423.  
  3424.                15. Repeat Steps 5 through Step 14.
  3425.  
  3426.                Optional computations:
  3427.  
  3428.                Compute geocentric and geodetic latitudes:
  3429.  
  3430.                    φ' = arctan(tan φ'),    φ = arctan(tan φ)
  3431.  
  3432.                Compute right ascension:
  3433.  
  3434.                    cos α = x/√(x2 + y2),    sin α = y/√(x2 + y2)
  3435.  
  3436.                    α = arctan(sin α /cos α)
  3437.  
  3438.                It is worthy  to note that  no trigonometric functions are
  3439.                actually  computed  with  this  algorithm,  yet  it  makes
  3440.                available the  sine and  cosine of  the geodetic latitude.
  3441.                Thus  the  rotation matrix,  required  for  computing wind
  3442.                velocity components, may readily  be constructed.  This is
  3443.                a great savings in  computer resources because this matrix
  3444.                must be computed for each integration step.
  3445.  
  3446.                The algorithm was tested  with latitudes ranging from -90°
  3447.                to 90°, and altitude ranging from -300 to 5 X 1012 feet at
  3448.                each latitude.  Miss distances of 10-1, 10-4 and 10-5 feet
  3449.                were used as  convergence criteria.  In  no case were more
  3450.                than three iterations required, and  no more than two were
  3451.                required for miss distances of 10-4 feet.
  3452.  
  3453.  
  3454.  
  3455.  
  3456.                Trajectory Maker                        Coordinate Systems
  3457.                Algorithm Supplement                                  6-15
  3458.                ──────────────────────────────────────────────────────────
  3459.                
  3460.  
  3461.  
  3462.                                        References
  3463.  
  3464.                1.  Anon., United States Navel Observatory and Royal
  3465.                    Greenwich Observatory, Explanatory Supplement to the
  3466.                    American Ephemeris and Nautical Almanac, Her Majesty's
  3467.                    Stationary Office, London, 1961.
  3468.  
  3469.                2.  O'Conner, J. J., Transformations Applicable to Missile
  3470.                    and Satellite Trajectory Computations, RCA
  3471.                    International Service Corporation, Technical Report
  3472.                    AFETR-TR-75_29, Patrick Air Force Base, Florida, July
  3473.                    1975.
  3474.  
  3475.                3.  Roy,A. E.,  The Foundations of Astrodynamics, The
  3476.                    Macmillan Company, New York, 1965.
  3477.  
  3478.                4.  Wertz, J. R., ed., Spacecraft Attitude Determination
  3479.                    and Control, D. Reidel Publishing Company, Boston,
  3480.                    1985.
  3481.  
  3482.  
  3483.  
  3484.  
  3485.                Trajectory Maker                              Bibliography
  3486.                Algorithms Supplement                                  7-1
  3487.                ──────────────────────────────────────────────────────────
  3488.                
  3489.  
  3490.  
  3491.                7.  BIBLIOGRAPHY
  3492.  
  3493.                Abell, G.,   Exploration of  the Universe,  Holt, Reinhart
  3494.                and Winston, New York, 1964.
  3495.  
  3496.                Anon., Space Flight  Handbooks, Volume I,  NASA SP-33 Part
  3497.                1, Office  of Scientific and  Technical Information, NASA,
  3498.                Washington D. C., 1963.
  3499.  
  3500.                Anon., United States Navel Observatory and Royal Greenwich
  3501.                observatory,  Explanatory   Supplement  to   the  American
  3502.                Ephemeris and  Nautical Almanac,  Her Majesty's Stationary
  3503.                Office, London, 1961.
  3504.  
  3505.                Anon., Orbital  Flight Handbook, Office  of Scientific and
  3506.                Technical  Information,  National  Aeronautics  and  Space
  3507.                Administration, Washington D.C., 1963.
  3508.  
  3509.                Bates,  R.,  and  D. Mueller,  J.  White,  Fundamentals of
  3510.                Astrodynamics, Dover Publications, Inc., New York, 1971.
  3511.  
  3512.                Battin,  B.,   An  Introduction  to  the  Mathematics  and
  3513.                Methods   of   Astrodynamics,    American   Institute   of
  3514.                Aeronautics and Astronautics, Inc., New York, 1987.
  3515.  
  3516.                Deetz,  C.  and  O. Adams,  Elements  of  Map Projections,
  3517.                Greenwood Press, New York, 1969.
  3518.  
  3519.                Felberg,  F., New  One-Step  Integration Methods  of High-
  3520.                -Order  Accuracy  Applied to  Some  Problems  in Celestial
  3521.                Mechanics, NASA Technical Report  NASA TR R-248, George C.
  3522.                Marshall Space Flight Center, Huntsville, Alabama, October
  3523.                1966, pp. 32-33.
  3524.  
  3525.                Gaposchkin, F.  M., ed.,  1973 Smithsonian  Standard Earth
  3526.                (III),  Smithsonian   Astrophysical  Observatory,  Special
  3527.                Report 353, Cambridge, Massachusetts, November 28, 1973.
  3528.  
  3529.                Hieb, H.,  Western Test Range  Geodetic Coordinate Manual,
  3530.                Part 1, Federal Electric Corporation, Feb. 1980.
  3531.  
  3532.                Ivy  L.  and H.  Reynolds,  On the  Sensitivity  of Impact
  3533.                Prediction   to  Methods   of   Interpolating  Atmospheric
  3534.                Density,  ITT Technical  Note, TN-76-1463,  Vandenberg Air
  3535.                Force Base, 1976.
  3536.  
  3537.                Maling,  D.,   Coordinate  Systems  and  Map  Projections,
  3538.                George Phillip and Son Limited, London, 1973.
  3539.  
  3540.  
  3541.  
  3542.  
  3543.                Trajectory Maker                              Bibliography
  3544.                Algorithms Supplement                                  7-2
  3545.                ──────────────────────────────────────────────────────────
  3546.                
  3547.  
  3548.  
  3549.                Mc  Cuskey,  S.,   Introduction  to  Celestial  Mechanics,
  3550.                Addison Wesley Publishing Company, Inc., Palo Alto,  1963.
  3551.  
  3552.                Menzner,  et.  al.,  Defining  Constants,  Equations,  and
  3553.                Abbreviated Tables of the  1975 U. S. Standard Atmosphere,
  3554.                National Aeronautics  and Space  Administration, Technical
  3555.                Report NASA TR R-459, Washington D. C., May 1976.
  3556.  
  3557.                O'Connor, J. J., Transformations Applicable to Missile and
  3558.                Satellite   Trajectory  Computations,   RCA  International
  3559.                Service  Corporation,   Technical  Report  AFETR-TR-75-29,
  3560.                Patrick Air Force Base, Florida, July 1975.
  3561.  
  3562.                Pearson,  F.,  Map  Projections: Theory  and Applications,
  3563.                CRC Press, Inc., Boca Raton, 1990.
  3564.  
  3565.                Ratcliffe, J., A.,  An Introduction to  the Ionosphere and
  3566.                Magnetosphere, Cambridge Univ. Press., 1972, p.31.
  3567.  
  3568.                Ridpath, I.,  Longman Illustrated  Dictionary of Astronomy
  3569.                and Astronautics, York Press, Essex, 1988.
  3570.  
  3571.                Roy, R.,  The Foundations of  Astrodynamics, The Macmillan
  3572.                Company, New York, 1965.
  3573.  
  3574.                Shanks,  B.  F., Solutions  of  Differential  Equations by
  3575.                evaluations of Functions, Math. Comp. 20 (1966) pp. 21-38.
  3576.  
  3577.                Steers,  D.,   An  Introduction   to  the   Study  of  Map
  3578.                Projections,  University  of  London  Press  Ltd., London,
  3579.                1962.
  3580.  
  3581.                Struick, D., Lectures  on Classical Differential Geometry,
  3582.                Addison-Wesley Publishing Company, Inc., Reading, 1961.
  3583.  
  3584.                Thompson,  R.  F.,  Spline   Interpolation  on  a  Digital
  3585.                Computer,  Goddard Space  Flight Center,  Technical Report
  3586.                X-692-70-261, Greenbelt, Md., July 1970.
  3587.  
  3588.                Waison, G.  R., Prediction of  Orbital Decay  Rates Due to
  3589.                Atmospheric Drag, TRW  Program Technical Report 3150-6020-
  3590.                RO-000, January 7, 1966.
  3591.  
  3592.                Wertz, J.  R., ed., Spacecraft  Attitude Determination and
  3593.                Control, D. Reidel Publishing Company, Boston, 1985.
  3594.  
  3595.                Whittaker, F. and G. Watson,  A Course in Modern Analysis,
  3596.                Cambridge University Press, Fourth Ed. 1963.
  3597.  
  3598.  
  3599.  
  3600.  
  3601.  
  3602.  
  3603.  
  3604.  
  3605.  
  3606.  
  3607.